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CALCULATION OF WATER DROP TRAJECTORIES TO AND ABOUT ARBITRARY 
THREE-DIMENSIONAL LIFTING AND NONLIFTING BODIES IN POTENTIAL AIRFLOW 


by Hillyer G. Norment 
Atmospheric Science Associates 

SUMMARY 

Computer programs are described by which trajectories of water drops 
can be calculated to and about three-dimensional bodies of arbitrary shape, 
which can have lifting surfaces, nonlifting surfaces or combinations of 
lifting and nonlifting surfaces. External potential airflow about a body 
is computed at subsonic air speed for any atmospheric conditions. Compressi- 
bility effects can be accounted for approximately. Flow into an inlet can 
be accommodated, provided that the intake flow rate is specified. To cal- 
culate water drop trajectories, experimentally derived relations between 
Reynolds and Davies numbers for water drops of all sizes, from the smallest 
cloud droplets to large raindrops, are used to represent effects of aero- 
dynamic drag on the particles during integration of the water drop equations 
of motion, and effects of gravity settling are included. A variable time 
step numerical integration method is used. 

The surface of the three-dimensional body is approximated by plane 
quadrilateral panels, over each of which a uniform potential source is 
assumed to be distributed. Source densities, lift vorticity and the re- 
sulting potential flow field are calculated by the Hess method. The codes 
are open-ended in their capacities to accommodate numbers of quadrilateral 
panels and other geometric complexities needed to define a body surface. 

The following seven codes are described: 

1, A code used to debug and plot body surface data. 

2. A modified version of the Hess lifting code which processes the 
body data and yields data required to compute flow velocities at 
arbitrary points in space. 


3. A code that computes flow velocities at arrays of points in three- 
dimensional space. 

4. A code that computes trajectories of water drops toward the body 
from arrays of initial points in space. 

5. A code that computes water drop trajectories and water drop 
fluxes to arbitrary target points. 

6o A code that computes water drop trajectories tangent to the body. 

7. A code that produces stereo pair plots that include both the 
body and trajectories. 

Code 1 and codes 4-7 are essentially the same as those described in NASA 
CR 3291. Code descriptions include operating instructions, card inputs and 
printouts for example problems. The FORTRAN IV source codes are available 
through COSMIC. 

Various tests of simulation accuracy are discussed, and, in general, 
accuracy is found to be acceptable. Calculated tangent trajectory results 
are compared with wind tunnel data and reasonable agreement is found. In 
comparison with the same experimental data, the new calculations show sub- 
stantial improvement over prior calculations. 
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INTRODUCTION 


A resurgence of interest in aircraft icing (ref. 1) has stimulated work 
on computer simulation of icing. Of critical importance to icing simulation 
are codes that compute trajectories of water drops, and other hydrometeors, 
to and about aircraft structures. We have recently published, in NASA CR 
3291 (ref. 2), a set of codes that compute water drop trajectories about 
arbitrary three-dimensional nonlifting bodies; the codes described here are 
similar to, in many cases are identical to, those described there. The new 
codes do everything the old ones do, plus account for effects on the flow of 
lift vorticity, if any. Another major difference is that the new codes use 
variable array dimensioning so that the storage requirements of the code can 
be tailored to the problem at hand. 

We distinguish two major categories of codes: flow codes and trajectory 

codes. The flow codes process data that describe the three-dimensional body 
and compute the fluid flow field around that body;* they are modified versions 
of the lifting flow codes developed by Hess (ref. 3). The trajectory codes 
use the flow codes' calculation results to compute trajectories of particles 
to and about the body. For aircraft icing studies, the body is, of course, 
an aircraft, the fluid is air and the particles usually are water drops. 

Table 1 identifies and briefly describes the executive codes in the two 
categories, and Table 2 does the same for the subroutine and function codes. 

As mentioned above, many of the codes are essentially the same as those 
described in NASA CR 3291 (ref. 2), and these are marked with asterisks in 
Tables 1 and 2. The user is cautioned to realize that while some of the new 
code subroutines may have the same names as their functional counterparts in 
the nonlifting code, it does not follow that they are the same. In fact they 
should be assumed to be entirely different unless they are marked with 
asterisks. 


*It is immaterial whether we consider the fluid to be stationary and the body 
in motion, or vice versa, but it is expedient here to consider the body 
stationary and the fluid in motion. 


3 


TABLE 1 


Code 

PBOXC* 

DUGLFT 

FLOPNT* 

Code 

ARYTRJ* 

CONFAC* 

TANTRA* 

STEREO* 

*This code 
described 


EXECUTIVE CODES 
A. FLOW CODES 

Description 

Processes and plots data which define the three-dimensional 
body. Used to debug and plot the body data. 

Processes three-dimensional lifting and/or nonlifting body 
data and prepares and stores data to be used by SR FLOVEL to 
calculate flow velocities as needed during trajectory 
calculations. 

Computes and prints flow velocities at user-specified arrays 
of points in space. 

B. TRAJECTORY CODES 

Description 

Computes trajectories, which begin at user-specified arrays 
of points in space, to and/or about the body. 

Computes trajectories from the free stream to user-specified 
points in space. Also computes particle concentration factors 
at user-specified points in space. (Concentration factor is 
ratio of particle flux at the target point to free stream 
particle flux. ) 

Computes trajectories tangent to the body which are initiated 
along user-specified lines in the free stream. (Tangent 
trajectories are those trajectories that barely miss inter- 
section with the body.) 

Prepares stereo-pair plots of the body along with particle 
trajectories. 


is essentially the same as the code of the same name that is 
in NASA CR 3291 (ref. 2). 
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TABLE 2 

SUBROUTINE AND FUNCTION CODES 
A. FLOW CODES 


Code 

Called [ 

AIJMX 

SIGMAL 

BVORTX 

CONTRL 

CKARRY 

CONTRL 

COLSOL 

SIGMAL 

CONTRL 

DUGLFT 

DATPRS 

INPUTL 

DKEKFK 

PISWIS 

FKUTTA 

BVORTX 

FL0VEL + 

TRAJEC* 

CONFAC* 

ARYTRJ* 

FLOPNT* 

HEADER 

CONTRL 

INPUTL 

NOLIFT 

LIFT 

VIJMX 

PNTVIJ 

DKEKFK 

UNIFLO 

SIGMAL 

AIJMX 

N I KMX 

VELOCY 

PRINTL 

INPUTL ^ 

CONTRL 

LIFT 

INPUTL 


— Description 

Computes the matrix, A _ , of dot products of source induced velocities with 
normal vectors to the on-body elements (ref. 3, eq. (7.12.1)). 

Calls PKUTTA or FKUTTA to calculate vortex strength per unit path length 
around the kth lifting strip, B^, for all lifting strips. 

Cross checks storage array capacities. 

Solves the linear equation matrices (ref. 3, eq. (7.12.5)) for element source 
densities. 

Controls flow of the modified Hess code which processes body surface data for 
use by the flow and trajectory codes. 

Translates, scales and rotates about the y axis surface description data 
immediately after input. 

Calculates Dj. £. and for use in calculating piecewise linear spanwise 
variation of B^ k '. (ref. 3, eq. (7.11.5)). 

Computes vortex strength, B (k) , for each lifting strip by the flow tangency 
method. 

Returns flow velocity for a given point in space. 


Writes a printout header. 


Inputs surface quadrilateral corner coordinates and controls computation of 
the geometric properties of the quadrilaterals. Also prints the first major 
DUGLFT output. 

Computes geometric properties of lifting quadrilateral elements (ref. 3, sec. 
7.2). 
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TABLE 2, cont. 


Code 
M I SI 

NEAR 

NEARF^ 

NIKMX 

NOLIFT 

PATPRS* 

PEADER* 

PICTURE 

PINPUT* 

PISWIS 

PKUTTA 

PNTVIJ 

PRINTL 

PSONST 

READ1 


Called By 

PKUTTA 

FKUTTA 

VFMLFT 
VFL I FT f 
SIGMAL 

INPUTL 

PINPUT* 

PINPUT* 

PBOXC* 

PBOXC* 

VMATRX 

BVORTX 

VMATRX 

VELOCY 

PISWIS 

SIGMAL 

COLSOL 

PKUTTA 

FKUTTA 

SUMSIG 

VELOCY 


Description 

(k) 

Linear equation solver. Used in calculation of vortex strengths, B v , of 
lifting strips. 


Computes source and vortex induced velocities from an element in a lifting 
strip using the near-field equations. 


Computes source and vortex induced velocities from an element in a lifting 
strip using the near-field equations. 

Computes the right hand sides of the A^ matrix, and n(°°\ which are the 

dot products of the onset flows with the unit normal vectors to the on-body 
elements, (ref. 3, eqs. (7.12.5)). 


Calculates geometric quantities for quadrilateral elements in a nonlifting 
section, (ref. 4, sec. 9.2). 


Translates, scales and rotates about the y axis surface description data 
immediately after input. 


Writes a printout header. 

Plots the body surface data. 

Processes input body coordinate data into quadrilaterals. Produces the non- 
lifting code "first output" (ref. 4, sec. 9.4). 

Calculates V^ and V^ according to equation (7.11.2) of reference 3, and 
calls DKEKFK and PSONST to calculate vortex induced onset flows. Used only 
for piecewise linear spanwise variation of vortex strength. 

Computes vortex strength, B^, for each lifting strip by the pressure equality 
method. 


If so requested 
velocities, V^ 


(for debugging purposes only), prints all 
and all vortex induced velocities, 


source induced 
and 


Prints the final output of the DUGLFT computations. 


Computes vortex induced onset flows when the piecewise-1 inear method of 
spanwise variation of vortex strength is used. 


Reads one singly subscripted array from a peripheral storage unit. 
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TABLE 2, cont. 


Code 

Called By 

Description 

READ3 

PNTVIJ 

STEPFN 

PISWIS 

PSONST 

UNIFLO 

AIJMX 

NIKMX 

PKUTTA 

VELOCY 

Reads three singly subscripted arrays from a peripheral storage unit. 

SETFLO 1 * 

FLOPNT* 

ARYTRJ* 

CONFAC* 

TANTRA* 

Reads DUGLFT output data stored on unit 18 that is required by SR FLOVEL 
for velocity calculations. If flow velocities are calculated by other than 
the Hess method, this code must be replaced with a dummy. 

SIGMAL 

CONTRL 

Controls calculation of the element source densities, and (ref. 3, 

eq. (7.12.5)). 

STEPFN 

VMATRX 

Computes vortex induced onset flows when the step function method of spanwise 
variation of vortex strength is used. 

ST0R18 

CKARRY 

Store control and quadrilateral geometrical property data on storage unit 18 
for use by the flow and trajectory codes, (see SR SETFLO.) 

SUMSIG 

CONTRL 

Computes the combined element source strengths, Oy (ref. 3, eq. (7.13.1)). 

UNIFLO 

VMATRX 

Stores uniform onset flow velocities for use in calculating element source 
densities. 

VELOCY 

CONTRL 

Computes the final velocity at the centroid of each element, and controls the 
final printout of the DUGLFT calculation, (ref. 3, eqs. (7.13.2) and (7.13.3) 

VFLIFT^ 

FLOVEL + 

Controls computation of velocities induced at a point in space by elements of 
unit source density and unit vortex strength in a lifting section. 

VFMLFT 

VIJMX 

Controls computation of velocities induced at all control points by elements 
of unit source density and unit vortex strength in a lifting section. 

VFMNLF 

VIJMX 

Computes velocities induced at all control points by elements of unit source 
density in a nonlifting section. 

VFNLFT + 

FLOVEL + 

Computes velocities induced at a point in space by elements of unit source 
density in a nonlifting section. 

VIJMX 

VMATRX 

Controls commutation of source induced velocities, V^, and vortex induced 
velocities, vj^ and vj?\ at the centroids of all elements. 

VMATRX 

CONTRL 

Subexecutive code for computation of induced and onset flow velocities at the 
centroids of all elements. 

WNEAR 

VFMLFT 

Computes vortex induced velocities from a wake element in a lifting strip 
using the near-field equations. 
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TABLE 2, cont. 


Code 

Called By 

Description 

WNEARF f 

VFNLFT + 

Computes vortex induced velocities from a wake element in a lifting strip using 
the near-field equations. 

WRITE! 

AIJMX 

NIKMX 

COLSOL 

PKUTTA 

FKUTTA 

SUMSIG 

Writes one singly subscripted array on to a peripheral storage unit. 

WRITE3 

VFMNLF 

VFMLFT 

STEPFN 

PISWIS 

PSONST 

UNIFLO 

Writes three singly subscripted arrays on to a peripheral storage unit. 
B. TRAJECTORY CODES 

Code 

Called By 

Description 

CDRR* 

PRFUN* 

PARTCL* 

Given Reynolds number, returns Davies number for a sphere. Used for water 
drops for which Reynolds number is less than or equal to 81.23. 

DVDQ* 

TRAJEC* 

Integrates particle equations of motion for each time step (ref. 7). 

FALWAT* 

PARTCL* 

Returns still -air, terminal settling speed for a water drop. Uses equations 
of Beard (ref. 13) . 

IMPACT* 

TRAJEC* 

Used in runs under control of CONFAC to adjust trajectory initial y,z 
coordinates to avoid impact on the body on the next trajectory after im- 
paction has occurred. This is a problem-specific subroutine that must be 
programmed by the user. (See p. 75.) 

MAP* 

CONFAC * 

Controls the iterative calculation of trajectories to a specified target 
point. 

MATINV* 

MAP* 

Linear equation solver. 

PARTCL* 

ARYTRJ * 
CONFAC * 
TANTRA * 

Read particle specification data and returns still -air, terminal particle 
settling speed and other particle data as required for the particular 
type of particle. This is a particle type-specific code. The version 
provided here is for water drops. 

POLYGO* 

CONFAC* 

Calculates area of a plane polygon of N vertices. Provides cross-sectional 
areas of particle flux tubes which are used to compute concentration factors, 
concentration ratios and collection efficiencies. 

PRFUN* 

TRAJEC* 

Given the particle Reynolds number, returns the factor which when multiplied 
by Vp - v a yields the first term on the right side of eq. (1). This is a 
particle type-specific function. The version provided here is for water drops. 
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TABLE 2, cont. 


Code 

Called By 

Description 

STRPNT* 

TANTRA* 

Specifies a curve in three-dimensional space on which lie the initial points 


of all trajectories used in computing a tangent trajectory to the body. Also 
specifies coarse and fine step sizes to be used in traversing the curve in 
search of the tangent trajectory, and it steps along the curve to define new 
initial trajectory points under control of TANTRA. The version supplied here 
uses straight line curves. 


T RAJ EC* 

ARYTRJ* 

TANTRA* 

MAP* 

Computes particle trajectories. (See p. 75.) 


TRANSF* 

PRFUN* 

MAP* 

Transforms coordinate system from the "flow system" 
or reverse. (See pp. 82-83.) 

1 to the "flux tube system". 

WCDRR* 

PRFUN* 

PARTCL* 

Given Reynolds number, returns Davies number for a 
where the Reynolds number is greater than 81.23. 

water drop. Used for case 


*This code is essentially the same as the code of the same name that is described in NASA CR 3291 (ref. 2). 
+Flow codes associated with SR FLOVEL rather than DUGLFT are so indicated by use of a superscript +. 
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METHODOLOGY 


THREE-DIMENSIONAL FLOW 

The methods and codes of Hess (ref. 3) and Hess and Smith (refs. 4,5) 
are used for calculation of lifting and nonlifting potential flow about 
arbitrary three-dimensional bodies. Lifting bodies (i.e., airfoils) alone, 
nonlifting bodies alone, or combinations of lifting bodies with nonlifting 
bodies (e.g., combinations of airfoils and fuselages) can be treated. Effects 
of flow into an inlet, for example an instrument aperature, can be accounted 
for provided the intake flow rate, in terms of fraction of free stream air 
speed, is specified. The method is restricted to subsonic airspeeds, but 
for free stream Mach numbers greater than 0.5, the Prandtl -Glauert method 
is used to correct approximately for compressibility effects. Since potential 
flow is computed, neither viscous effects nor turbulence are treated. 

The code requires input of a digital description of the body surface, 
and for purposes of organizing the data as well as for computing flow, the 
body surface is partitioned into sections which are designated as either 
lifting or nonlifting. In either case, the surface is represented by con- 
tiguous, plane quadrilateral panels (usually called elements). (See Fig. 1.) 
For nonlifting sections there are few restrictions on the manner in which 
the elements can be arranged to represent the surface other than those re- 
quired for organization. Lifting sections are restricted as follows: each 

must consist of strips of elements, the strips being oriented parallel to 
the chordwise direction of the airfoil; each strip must have the same num- 
ber of elements; and wake elements must be included after the trailing 
edge of each strip. Details are given below on pp. 20-27. Both lifting 
and nonlifting portions of the body may be described by more than one section. 

Each on-body element (which is in the flow) is taken to be a potential 
flow source. The source is a distributed one, with the distribution being 
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a. Forward fuselage of the C130E. 


Leading Edge 




uniform over the surface of the element, and each element, for example the 
jth s -j s characterized by a unique source density, a-. In addition, each 

J 

strip of elements in a lifting section is characterized by having a unique 
value of lift vorticity associated with it. This quantity, for example for 
the lifting strip, B^, represents vortex strength per unit path length 
around the strip (see Fig. 4 and its discussion), and it represents the sum 
of contributions from all panels in the strip. Velocities induced by these 
vorticities are treated as onset flows. Thus, there is an onset flow from 
each lifting strip plus the free stream onset flow. It is necessary to com- 
pute an independent source density for each of these onset flows for each 
on-body quadrilateral panel: if there are N on-body panels, K lifting strips 

and one free stream flow, N ( K + 1 ) values of o must be computed. Source den- 
sities are determined by solving large systems of linear equations that 
represent the effects of all onset flows on all panels, plus the mutual inter- 
actions of all distributed sources, under boundary conditions of zero flux 
through the centroid (also called control point) of each on-body panel, or 
specified fraction of free stream flux through each inlet panel.* 

Determination of vortex strengths requires an additional constraint, the 
Kutta condition, and this is supplied by user-selection of one of two optional 
methods which are designated as "flow tangency" and "pressure equality." 
Application of these options is discussed below on pp. 38-39. 

Lift vorticity is computed by a novel method developed by Hess (ref. 3). 
To circumvent problems that have been found to result from use of vortex 
filaments in prior work, and to ensure that potential flow results from the 
vorticity distribution and that individual infinitesimal vortex lines either 
form closed curves or go to infinity, Hess has developed a method by which 
vortex sheets on the body and wake surfaces can be expressed in terms of 
dipole sheets on the same surfaces. Hess summarizes the method as follows: 

"A variable-strength dipole sheet is equivalent to the sum of: 

(1) a variable-strength vortex sheet on the same surface as the 
dipole sheet whose vorticity has a direction at right angles to 


*An inlet orifice is paneled just as is an impermeable body surface (see Fig. 6). 
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the gradient of the dipole strength and a magnitude equal to the 
magnitude of this gradient, and (2) a concentrated vortex filament 
around the edge of the sheet whose strength is everywhere equal to 
the local edge value of dipole strength." 

Mathematical details are given in Appendix A of reference 3. 

For particular body geometry and orientation relative to the free 
stream, the source densities and vortex strengths are calculated only 
once, and then these can be used to calculate flow velocity at any space 
point exterior to the body. The primary functions of the DUGLFT codes are 
to calculate the ck and and store these quantities, along with other 
requisite data, for use by subroutine FLOVEL in calculating flow velocities. 
Subroutine FLOVEL is called as needed by programs TRAJEC, CONFAC, ARYTRJ 
and FLOPNT to provide flow velocities for trajectory and flow velocity array 
calculations. 

In calculating each flow velocity, contributions from all quadrilateral 
elements are summed. There are three sets of algorithms for computing con- 
tributions from individual elements: (1) for elements that are close to the 

calculation point, detailed calculations are used that account for exact 
element geometries, (2) for elements at intermediate distances multi pole 
expansions are used, and (3) for remote elements the point source approximation 
is used. Mathematical details are given in references 3-5, with emphasis on 
lifting flow in reference 3 and emphasis on nonlifting flow in references 4 
and 5. The reader is strongly urged to study these references closely before 
attempting to use this code. Reference 6 consists of a code users manual 
for the lifting flow calculations described in reference 3. 

Calculation accuracy is discussed below in the section on Validation 
(p. 95). Of course accuracy also depends on the fineness of resolution of 
the element description of the body, and naturally some compromise is called 
for. The smaller the elements the finer the resolution, and the fewer of 
them for which the most exacting of the three algorithms must be used. On 
the other hand, the number of elements increases inversely as the square of 
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their linear size. In past studies on airplanes we have used the following 
paneling criteria: For those parts of the airplane traversed by particle 
trajectories, we try to keep the element edges between 6" and 8" in length. 
Where allowed by simplicity of surface shape, remote elements can be larger. 
Remote downstream complexities of shape are ignored or treated approximately. 
For example, if interest is confined to the forward fuselage, then the re- 
mainder of the fuselage can be represented as a cylinder of constant cross- 
section which is extended to approximately five times the length of the nose 
section (as recommended by Hess and Smith, ref. 5), and the wings can be 
ignored entirely. 

The following are basic requirements of the method that apply to all 
calculations : 

1. A uniform, unit-speed free stream approximately in the direction 
of the positive x axis. 

2. Normalization of all velocities to be consistent with the unit 
free stream speed. 

3. Normalization of all distances by a user-specified characteristic 
dimension of the body. 

Surface point coordinates may be recorded in any convenient units and can be 
appropriately translated and scaled, to meet requirement 3 above, during 
processing via use of SR's PATPRS and DATPRS. These subroutines also allow 
rotation of the body about the y axis to adjust attitude angle.* The 
coordinate system used for the calculations is described on pp. 19-20. 

The unit free stream speed is assumed by program DUGLFT , and the distance 
normalization, if required, is done during preliminary data processing as 
indicated above. For trajectory calculations, the user specifies the true 
free stream speed and the normalization length, and the codes automatically 
handle any additional normalizing or scaling that is required. 


translating, scaling and rotating are controlled via parameters IPROS and 
data cards 4 and 5 of programs PBOXC and DUGLFT, respectively. 
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PARTICLE TRAJECTORY CALCULATION 


We assume that the bulk air flow is not perturbed by the particles. 
Moreover, since particle density is large compared to that of air, we can 
neglect buoyancy and inertial reaction of the fluid to obtain the three- 
dimensional, normalized equation 



Non-dimensional quantities are: 


t 

T 

n d = c d n r 

N f = V 2 /(Lg) 

N r = fii |v - v | V 
R n a p 1 


particle and air velocities 

still-air, terminal settling speed of the particle 

Unit vector in the z (upward) direction 

time 

Davies number 
Froude number 
Reynolds number 
Particle drag coefficient 


Dimensioned quantities are: 

6 particle diameter 

p air density 

n air viscosity 

g gravity acceleration constant 

V free stream airspeed 

L a characteristic dimension of the body 
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Here length is normalized by dividing by L, velocity by V and time by L/V. 

N n and N D are for still-air, terminal particle settling. 

We initiate the calculation far enough upstream to be essentially 
beyond the influence of the body where we can take v = v - l<v . We 
compute N r from these data, calculate Np from N R using the relations 
discussed in the next section, and proceed straightforwardly with a 
numerical integration of eq. (1). The integration is done via use of 
the code DVDQ of Krogh (ref. 7). This code uses an Adams-type predictor - 
corrector algorithm with variable time step. It also tests for compu- 
tational stability and loss of accuracy via roundoff error. It was tested 
by Hull, et al (ref. 8), along with a number of other codes and found to 
be most efficient in terms of numbers of function evaluations (flow 
velocities) required. 

For a particular case, time required for trajectory calculation is 
largely dependent on the number of elements and the number of velocities 
required. On the CDC 6600 computer, one nonlifting velocity calculation 
requires on the order of 0.15 second for a typical problem. The number of 
velocities required per trajectory varies from about 60 to 300. A typical 
number of trajectories required is 25. Thus, computing time, even on a 
large computer, can be considerable. Computing times required for the test 
problem and validation runs are given below on pp. 103-106. 

AERODYNAMIC DRAG ON WATER DROPS 

Davies (ref. 9) shows that still -air terminal settling of spheres can 
be generalized in terms of the dimensionless numbers N R s and N^ s - Over 
the range from the smallest spheres, which settle under viscous flow con- 
ditions and obey Stokes law, to spheres much larger than of interest here, 
and for any Newtonian fluid, a reproducible single-valued relationship 
between N D _ and N n exists. Furthermore, N n is independent of settling 
speed, being a function of fluid and sphere properties only; thus for given 
sphere and fluid, N D and hence V can be calculated. Polynomials by which 
N d can be computed as a function of N n were derived by Davies from a 

K j S U 5 5 

composite of many sets of experimental data. 
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Since the work of Davies it has been found repeatedly that this treat- 
ment is applicable to particles of other shapes, provided settling is 
steady and particle orientation is stable. 

For the trajectory calculations required here, the problem must be 
turned around. In addition to gravity settling, there is a particle vel- 
ocity component (relative to air) caused by the disturbance of the passing 
airplane. At any time step in the numerical integration of equation (1), 

v - v (and hence N n ) is known, and N n must be determined. For viscous 
a p R D 

motion (i.e., Stokes flow, where N R < 1) Np = 24 N R and equation (1) can be 
integrated without question. However, for larger N R the steady-state drag 
data determined experimentally for terminal settling must be used to com- 
pute accelerative particle motion. 

Experimental measurements by Keim (ref. 10) and a theoretical analysis 
by Crowe, et al. (ref. 11) indicate that if the acceleration modulus, 


n a 




5 


is smaller than about 10 -2 , steady-state drag coefficients can be used 
without significant error to compute accelerative motion. N^ rarely 
exceeds 10“ 2 in our trajectory calculations. 

For water drops small enough to be essentially spherical (N R < 81.23) 
we calculate Np from a polynomial function in N R derived from Davies data 
(ref. 9). (Function CDRR) For larger drops (N R > 81.23), which have a 
flattened, non-spherical shape, we calculate Np from polynomials in N R 
derived from the water drop data of Gunn and Kinser (ref. 12). (Function 
WCDRR). 

Still-air, terminal settling speeds for water drops are computed via 
use of Beard's equations (ref. 13). (SR FALWAT) 

Water drops of any size, from submicron to the breakup size at about 
8000 ym diameter, can be handled by these methods. However, the user 
should be aware that computation time goes up as droplet diameter goes 
down, and the time required for drops of diameter 1 ym or less may be large. 
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FLOW CODE DESCRIPTIONS 


PROGRAM PBOXC 
General Discussion 

This program is derived from the Douglas Aircraft Company code BOXC 
which was developed by Hess and Smith (ref. 4) and is identical to the code 
described in reference 2. It processes and produces CALCOMP plots of the 
three-dimensional body surface description data and is used primarily to 
debug these data. Processing and printing go as far as the "first output" 

(ref. 4, sec. 9.4). A secondary use is to store the body surface data such 
that it can be retrieved later and used by PGM STEREO to plot the body along 
with trajectories stored by one of the trajectory codes. 

The surface of a general three-dimensional body is defined in terms of 
"rows" and "columns", the so-called m and n-l’nes, of coordinates of points 
on the surface as described below. The m and n-lines of points are combined 
by the code to form quadrilateral elements, or panels, such that when con- 
sidered together they represent a reasonable approximation to the surface 
(Fig. 1). Adjacent panels should be contiguous, or as nearly contiguous as 
possible. The data for general bodies may be scaled and translated in the 
three coordinate directions, and rotated about the y axis (see the next 
section) prior to processing. 

This code also has the capability to generate ellipsoids of prolate, 
oblate or general shape, with the only restriction being that their major 
and minor axes lie on the coordinate axes,* 

When the user elects to prepare plots of the body, the code automatically 
prepares a number of plots, each from a unique viewing angle, the number vary- 
ing according to symmetry. For an asymmetric body fourteen plots are pre- 
pared. These consist of the six views from both directions along each 
coordinate axis, and the eight plots from 45 degree angles in each octant. 


*PGM DUGLFT does not have this capability. 


T8 


For a body with one plane of symmetry nine plots are prepared, for two 
symmetry planes six plots, and for three planes four plots. The user is 
urged to make liberal use of the plots to find errors in the body data. 

Coordinate System 

Figure 2 shows an example of a cartesian coordinate system that would 
be appropriate for description of the body surface and for calculation of 
flow velocities and particle trajectories. Critical directions are those 
of the free stream and the gravity vector. During particle trajectory 
calculations the z axis must be directed vertically upward such that gravity 
settling of particles is along the negative z axis. The x axis typically is 
taken to be in the plane that contains the z axis and free stream vector, 
with its positive direction the same as that of the free stream vector 
projection on the x axis. The y axis then has the direction required by a 
right-handed coordinate system. 



Figure 2. Coordinate system suitable for body surface description 
and particle trajectory calculations. 
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For a body such as an airplane, the orientation shown in Figure 2 is 
the proper one to use in setting up the surface digital description. Any 
convenient units and any origin can be used since these can be adjusted at 
data processing time by application of a coordinate translation and scaling 
capability (SR PATPRS and DATPRS). This capability also allows rotation 
about the y axis such as to adjust angle of attack. As discussed below on 

p. 37, angle of attack also can be specified arbitrarily during the DUGLFT 
calculations, but if particle trajectories are to be computed, the user must 
keep in mind that the negative z axis must have the direction of the gravity 
vector. 

As noted in the discussion of equation (1), all distances are normalized by 
dividing by a characteristic dimension of the body. Examples of appropriate 
characteristic body dimensions are fuselage diameter or airfoil chord length. 
Therefore, when the coordinate data are scaled by PATPRS or DATPRS, it is 
appropriate to define the scale factors to be the reciprocal of this charac- 
teristic length. 

If the body has one or more planes of symmetry, these can be used by 
the code as reflection planes to generate the complete body from the basic 
unit structure as described in the next section. However, the origin of 
the coordinate system must be adjusted (via coordinate translations) such 
that the symmetry planes are coincident with the y = 0 and z = 0 planes. 

For example, the symmetry plane of the airplane in Figure 2 must be coinci- 
dent with the y = 0 plane. More details are given in the next section. Again, 
the restriction that the z axis must be directed opposite to the gravity 
vector may complicate matters. For example, if a body normally has a 
horizontal (z = 0) symmetry plane, but for some particular case must be 
oriented with an angle-of-attack to a free stream that is horizontal rela- 
tive to the ground, the symmetry could not be used since the z axis would 
no longer be normal to the symmetry plane. 

Symmetry Planes 

The number of symmetry planes is specified via input parameter NSYM 
which has allowed values of 0,1, 2, 3. Only "plus" symmetry planes are 
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allowed. (See p. 37 for a description of plus and minus symmetry planes.) 
Symmetry planes are specified and used for reflection according to a fixed 
procedure: If there is one symmetry plane it is the y = 0 plane. If there 

is more than one symmetry plane (NSYM > 1) the first reflection is always 
across the y = 0 plane, followed by reflection across the z = 0 plane, which 
latter is always the second reflection plane. Specification of three sym- 
metry planes (the third symmetry plane is always the x = 0 plane) is used 
only to restrict the number of perspective view plots, and in generating 
ellipsoids for the case of I FLAG = 1 (pp. 27-29); surface point coordinates are 
not reflected across the third plane during calculation of the plots, and 
a third reflection plane is not allowed by the DUGLFT code. 

In summary, the reflection symmetry planes are exercised as follows: 

Symmetry 
Plane 

y = 0 

z = 0 

For example, if NSYM = 1, for each point with coordinates (x,y,z), another 
point with coordinates (x,-y,z) is created. If NSYM = 2, for each point 
with coordinates x,y,z, three additional points with coordinates (x,-y,z), 
(x,-y,-z) and (x,y,-z) are created. 

Only the primary data points should be input if symmetry is to be used. 
If reflected as well as primary data are input, and both are reflected by 
symmetry, the calculation results will be in error. 

Read the Uniform Onset Flow section, p. 37, before proceeding with 
use of symmetry planes. 




The user must examine the body, or drawings of it, and devise a layout 
plan for subdividing its surface into sections that are compatible with the 
requirements of m-line, n-line surface point input for lifting and nonlifting 
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surfaces while providing elements of appropriate size which cover the surface 
without leaving gaps or introducing unwanted discontinuities. Also a co- 
ordinate system must be established, but this can be manipulated at processing 
time by use of the scaling and translation capability of the code as noted above. 

The important thing is to understand the requirements of the m and n-line 
input. Here we give a brief summary of the requirements; the user is 
encouraged to carefully study sec. 9.1 of reference 4, and secs. 7.1 - 7.3.1 
of reference 3 to obtain a thorough understanding of them. 

Points which define the corners of the quadrilateral elements are labeled 
with integers m and n which identify hypothetical "rows" and "columns" on 
which they lie. The integers m and n are not input to the computer; they 
are used for data organization and sequencing only. 

To ensure a proper computation, the rows and columns must be organized 
by the following rule: If an observer is located in the flow and is oriented 
so that locally he sees points on the surface with m values increasi ng up- 
ward (or forward), he must also see n values increasing toward the right. 

The surface of any part of the body may be described by more than one 
section, each of which must be independent. That is, all quadrilaterals in 
each section must be closed. Where an edge of a section is contiguous with 
another, the input for each section must define the common edge, though 
they need not use the same points on the edge. 

Figure 3 illustrates a nonlifting surface description that is subdivided 
into four sections. Note how the sectioning can be used to change resolution 
or to deal with structural complexities. Fully general quadrilateral shapes 
are acceptable in nonlifting sections, including triangles, in which case 
two points may be coincident or one of the four points may lie on a triangle 
edge. 

Coordinates are punched into cards, one point per card; also in each 
card is punched the integer parameter STAT which is used to identify the 
m,n status of each point. All points in a section are ordered in the se- 
quence (m,n): 

(1,1), (2,1), (3,1),..., (1,2), (2,2), (3, 2),. ..(1,3), (2,3), (3,3),... 
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Figure 3. Plan view of the input points on a body divided 
into sections. (From ref. 4) 


The STAT parameters are as follows for each section, whether lifting or 
nonlifting: 

(m,n) STAT 

( 1 , 1 ) 2 

(l.rtfl) 1 

all others 0 or blank 

For the last card of the last section, STAT = 3. 

Input order of sections is immaterial, but within sections, the data 
must be ordered according to the underlined rule given above. 
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Surface Description Data for Lifting Section - General 


For lifting sections there are certain restrictions and requirements 
beyond those for general nonlifting sections. The n- lines in a lifting 
section must be approximately parallel with each other, each must have the 
same number of points, and each must lie roughly in the free stream direc- 
tion. The code assumes a trapazoidal shape for each lifting element, with 
the parallel sides of the trapazoids along the n-lines. (An n-line is a 
line along which n is constant; in this case it is chordwise.) In addition 
to the elements which define the lifting surface, elements also must be 
defined for the wake. These wake elements are constructed similar to the 
body surface elements, the first is contiguous with the body elements 
along the trailing edge, but they extend aftward of the body to represent 
an infinitesimally thick vortex sheet wake. Wake elements are assigned 
vorticity but not source density since, of course, they are not on the body 
surface. 

Figure 4 illustrates how an airfoil surface and wake should be parti- 
tioned into elements for a typical situation (i.e., a wing with trailing 
edge normal to the free stream flow or tapering backward, and a zero or 
nose-up angle of attack). In any situation n varies spanwise while m 
varies chordwise. For the typical situation, the m index begins at the 
trailing edge .proceeds forward along the lower surface, around the leading 
edge, aftward along the upper surface to the trailing edge and then into 
the wake. To satisfy the "underlined rule" on p. 22 for organizing input 
data, it is common practice to define the wing on the starboard side of the 
airplane with n = 1 at the tip and proceed inward toward the root. 

For cases where the trailing edge tapers forward at a significant 
angle, or negative lift is expected, it would seem best (ref. 3, secs. 

6.5 and 8.4) to reverse the above procedure. That is, begin at the trailing 
edge as before, but traverse the upper surface first, around the leading 
edge, back across the lower surface and finally into the wake. On a star- 
board wing, this reverse procedure requires that n values increase from 
root to tip, and the code allows this order of input as well as the more 
usual tip to root order. 
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DIRECTION 


Figure 4. Organization of m and n lines in a lifting section. A lifting 
strip is delineated by sequential n lines, and extends over the 
complete circuit from m = 1 at the trailing edge, along the 
underside to and around the leading edge, back to the trailing 
edge, and finally to the furthest aftward extent of the wake. 
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Wake elements may be defined downstream as far as desired. The user 
has the option, which is recommended, to specify, via election of the semi- 
infinite last wake element option (see p. 35), that the code extend the 
trailing edge of the last defined wake element to infinity. Thus, the user 
should specify enough wake elements to account for whatever curvature is 
considered necessary, plus one final element, and exercise the semi-infinite 
wake element option; two or three wake elements should be sufficient. 

The considerations above of wake structure and position are mitigated 
by the following. First it is not known a priori exactly how a wake for a 
particular case should be configured. Second, it has been found that the 
flow on a lifting body is insensitive to the position of its own wake, and 
thus, for ordinary moderate values of: angle of attack, trailing edge angle 

and degree of camber, any reasonable wake gives a satisfactory solution 
(ref. 3, sec. 8.5). 

When a wing is subdivided into several lifting sections, the organi- 
zation proceeds in the same directions within each section, and in the 
typical situation begins with the outermost section and proceeds inboard 
in orderly succession of sections. 

Hess has studied flow calculations around an isolated wing to determine 
the number of lifting strips and elements per strip that are needed for 
accurate results (ref. 3, sec. 8.1). He found that thirteen lifting strips 
per half wing and about thirty elements per strip, equally divided between 
the upper and lower surfaces, were satisfactory. 


Surface Description Data for Lifting Sections - Extra Strips And Ignored 
Elements ~ 

Two special situations arise from intersection of nonlifting structures 
with lifting sections. First there is the case where the end of a lifting 
section abuts a nonlifting structure along the entire length of both its 
upper and lower surfaces. An example is the intersection of a wing with a 
fuselage. The lifting section should not simply terminate at the inter- 
section because the local lift does not fall to zero there, and in addition, 
termination would result in a concentrated vortex filament on the surface. 
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To remedy this, an "extra strip" is added to the end of the lifting section 
which extends into the nonlifting body. This extra strip is geometrically 
similar to the others in the section and includes a wake. (Of course, 
element source densities are not computed for elements in extra strips.) 

In the case of intersection with a fuselage, the extra strip extends the 
lifting section to the airplane symmetry plane. These extra strips must be 
specified by the user in the same way that ordinary strips are specified. 

The second situation arises when a nonlifting structure covers part of 
a lifting section, usually over only a portion of a lifting strip, such as to 
not break the continuity of the trailing edge. Examples would be engine 
cowlings and wing pylons. The affected elements are designated as "ignored 
elements", but the vortex distribution continuity is maintained over them to 
avoid numerical problems. Thus as regards the vortex distribution, the 
lifting surface is taken to extend through the nonlifting structure as 
though the nonlifting structure did not exist. Of course, source densities 
are not computed for the ignored elements, and they contribute nothing to 
the nonlifting flow. 

Special extra strips are required at junctions of independent, con- 
tiguous lifting sections when the piecewise linear option is used to 
smooth the spanwise vortex distribution. However, these special extra 
strips are specified by the code on the basis of information supplied via 
the NLINE1 and NLINEN parameters (pp. 40-42), and the user does not need to 
otherwise concern himself with them. 


Surface Description Data for Generated Ellipsoids (I FLAG = 1 or 2) 

Ellipsoids are generated by specifying the semi-axis lengths 
B and C (A = 1 always), and by specifying the numbers of "latitudinal" 

and "longitudinal" (Fig. 5) element divisions, NLM1 and MM IN respec- 
tively. 
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S ^ — . ellipsoid cross section 

in the x = x x plane 
(n constant) 


9 is in the y = 0 plane 
<t> is in the x = xi plane 


m-lines run in the latitudinal direction from e = 0 to 0 = ■fi- 
n-lines run in the longitudinal direction from <j> = 0 to <J> = 2 


Figure 5. Definition of angles 0 and <|>, and m and n-line directions 
used by PBOXC for generation of ellipsoids. 


There are two modes for specification: 


Mode 1. I FLAG = 1, NSYM = 3 

All three symmetry planes are used and NLM1 and 
MMIN are specified for one octant only. Element 
increments are computed for NLM1 and MMIN equal 
increments in angles 0 and <|> (Fig. 5). 


Mode 2. I FLAG = 2, NSYM = 2 

Only two symmetry planes are used, and (x, z) 
values in the y = 0 plane must be input for 
-1 < x < 1, beginning at (1,0) and proceeding 
to (-1,0) for either all positive z or all 
negative z (i.e., for 180° in angle e). (The 
code automatically ensures that the "underlined" 
input rule is obeyed.) Thus, NLM1 must be speci- 
fied for the entire x axis, but MMIN is for one 
octant only as for the other option, and element 
increments in the "longitudinal 11 direction are 
created at equal increments of the angle <f>. 

Body surface data for generated ellipsoids cannot be plotted nor can 
the data be translated, scaled, and rotated by subroutine PATPRS. DUGLFT 
does not have a capability to generate ellipsoids. This capability is used 
only by the nonlifting code described in reference 2. 
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Printed Output 


The printed output is the result of the first stage of surface data 
processing (ref. 4, sec. 9.4) and is obtained by specifying IPRNT = .TRUE, 
in the input. For each quadrilateral element on the surface it consists of: 

1. Coordinates (X,Y,Z) of the four points on a quadrilateral in the 
order 



m 

4 


>. n 


around the quadrilateral. 

2. Components (NX,NY,NZ) of the unit normal vector to the plane of 
the quadrilateral. This vector should point toward the exterior 
of the body rather than toward its interior. If it points in the 
wrong direction the data have been input in violation of the 
"underlined rule" on p. 22 and the data must be reordered. 

3. Coordinates (NPX,NPY,NPZ) of the quadrilateral null point (ref. 4, 
sec. 9.3). 

4. The common projection distance (D) of the four input points into 
the plane of the quadrilateral. (The four points from which a 
plane quadrilateral is formed do not in general, and need not, lie 
exactly in a plane.) 

5. The maximum diagonal length (T) of the quadrilateral. 

6. The area (A) of the quadrilateral. 
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Additional output appears for certain abnormal quadrilaterals. If the 
integer 1 or 2 appears at the far right of the page, they indicate the 
following conditions: 

Integer 1. The null point was found to lie outside of the 
quadrilateral. The coordinates listed are for 
the quadrilateral centroid. 

Integer 2. The iterative procedure used to determine the null 

point did not converge and thus the null point is only 
approximate. 

(ref. 4, sec. 9.3). 


Subroutines Required 

PINPUT, PICTUR, PEADER, PATPRS, plotting subroutines. 

External Storage Units 

Units 1, 5 and 6 are the system punch, input and print units respectively. 
Unit 8 temporary storage. 

Unit 9 storage for surface data to be used later for plotting by 
PGM STEREO. 
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PBOXC Card Input 


Card 

No. Variables and Format 

1 HEDR( 15 ) , I FLAG, NSYM, 

KMACH, KASE, (15A4, II, 
10X, II, IX, II, 2X, A4) 


2 MACH, (F10.6) 


HEDR (Cols. 1-60) 
I FLAG (Col. 61) 

I FLAG = 0 

I FLAG = 1 


Description 

Hollerith run identification 

Body surface description input control 

Input data for a general, three-dimensional 
body (See pp. 21 ff.) 

Generate an ellipsoid using the mode 1 option, 
with three reflection planes. (See p. 29.) 

Be sure that NSYM = 3. 


I FLAG = 2 


NSYM (Col. 72) 
KMACH (Col. 74) 


Generate an ellipsoid using the mode 2 option, 
with two reflection planes, and input x,z 
coordinates for the ellipsoid via cards no. 5C. 
(See p. 29.) Be sure that NSYM = 

Note: The generated ellipsoid capability 

( IFLAG f 0) applies only to the non- 
lifting code (ref. 2). 

Number of data reflection planes. Limited to 
values 0,1, 2, 3. (See p. 20.) 

Always zero or blank. 


KASE (Cols. 77-80) Hollerith body identification. 


Mach number This card is input only if KMACH / 0 on card 1. 


3 IPROS, IPUNCH,IPRNT, 

I PICT , ICRT, (5L1) 


4. ANGLE, XSCALE , YSCALE, 

ZSCALE , XTRANS, YTRANS , 
ZTRANS , (7F10.0) 


Logical variables which cause the following if true: 

IPROS Body surface data for a general body are to be trans- 
lated, scaled and rotated about the y axis before 
processing, and card 4 is to be input. 

IPUNCH Body surface data are copied to the system punch unit. 

(After translating, scaling and rotating about the y axis.) 

IPRNT Print of the output described on p. 30. 

I PICT Body surface data for a general body are plotted. 

ICRT Plotting is via CRT. If ICRT is false, plotting is via 

pen and ink. 


This card is input only if IPROS (card 3) is true. 

ANGLE Angle (degrees) that the body is rotated about the y axis. 

A positive value causes a counterclockwise rotation from 
the aspect of a viewer looking down the positive y axis 
toward the origin. (Note: For a nose-up airplane angle 

of attack, ANGLE is negative.) 

XSCALE, Scale factors to be applied to surface point x, y and z 

YSCALE, coordinates respectively after translation. Default 

ZSCALE values are unity. 

XTRANS, Translations to be applied to surface point x, y and z 

YTRANS, coordinates before scaling. 

ZTRANS 


5A X, Y ,Z ,STAT, XX, YY , 

ZZ, STATT, 

(3F10.0 , 1 2/3 F10. 0, 12) 


Cards 5A apply to general bodies (IFLAG = 0, see pp. 21 ff). 

X,Y,Z and Are coordinates of points used to define the body 
XX, YY, ZZ surface. 


STAT Are point status integers. Allowed values are 0, 1, 

STATT 2, 3. The meaninqs of these values are: 

(Col. 32) 0 This point is on the same n line as the last point 

1 This point starts a new n line 

2 This point starts a new section 

3 This is the last point in the input. 
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PBOXC Card Input, cont. 


Card 

No. Variables and Formats Description 


Note: For the last coordinate card STAT or STATT = 3. A Mank 

card should follow this if there is an odd number of body 
surface point cards. 


5B NLM1, MMIN, B, C, 

(215, 2F10.5) 


Card 5B applies to generated ellipsoids (IFLAG > 0) 
NLM1 Number of "latitudinal" element divisions 
MMIN Number of "longitudinal" element divisions 
B y semi-axis of the ellipse 

C z semi -axis of the ellipse 

(See p. 27 Modes 1 and 2, and Fig. 5.) 


5C x ls zi, x 2 , z 2 , 

**' X NLMl+r Z NLM1+1 
(8F10.0) 


Cards 5C apply to qenerated ellipsoids for which the x,z coordinates 
are input (IFLAG = 2, NSYM = 2). 


(x-, z.) are coordinates in the y = 0 plane, beginning at (1,0) and 
prbceeding to (-1,0), that define the "latitudinal" element subdi- 
visions. (See p. 29 Mode 2, and Fig. 5.) 


LINE1, LINE2, (7A6/7A6) Cards 6 and 7 are read only if ICRT is true (card 3). These are 

two lines of 42 columns each of Hollerith labeling for a micro- 
fiche film. 
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DU6LFT 


General Discussion 

Program DUGLFT is an extensively modified version of the lifting code 
developed by Hess (ref. 3) and documented for users by Mack (ref. 6). Modi- 
fications were required to accomplish the following: 

1. Provide data output (on storage unit 18) for use by SR FLOVEL 
for calculation of flow velocities at arbitrary space points. 

2. Provide capability to treat nonlifting as well as lifting bodies, 
and combinations of lifting and nonlifting bodies. 

3. Provide capability to translate. scale and rotate body surface 
coordinates. 

4. Provide capability to treat flow inlets. 

5. Provide approximate correction for compressibility effects at 
high subsonic free stream speeds. 

6. Generalize array capacities in the code such as to make the code 
open ended with regard to its capacity to handle numbers of: 
quadrilateral elements, sections, lifting and nonlifting strips, 
etc. 

The discussions to follow are intended to provide guidance to code users. 
An understanding of the method, its strengths, weaknesses and range of poten- 
tial applications can only be obtained by thorough study of references 3, 4 
and 5. 

Body Surface and Make Description 

Coordinates of corner points of the plane quadrilateral elements which 
described the body surface and wake are organized and input via cards no. 12 
as described above for PGM PBOXC (pp. 19-27). 
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Semi-Infinite Last Wake Element Option 


By specifying LASWAK = .TRUE, on card 7, the user specifies that the 
semi-infinite last wake option be exercised. In this case, the trailing edge 
of the last wake element in each lifting strip is automatically extended to 
x = co by the program. Thus, the wake vorticity also extends aftward to 
infinity, as is consistent with theory, and the user is encouraged to exercise 
this option. 

All four corner points, including points that define the trailing edge, 
of each of the semi-infinite last wake elements must be specified by input 
(cards no. 12) just as is done for all other elements. This is done such 
that a set of finite geometric quantities, including a control (i.e., cen- 
troid) point, can be calculated. For example. Fig. lb shows three wake 
elements on each lifting strip of the C130 wing. If the semi-infinite last 
wake element option is exercised, geometric properties would be calculated 
for the last (most aftward) of each of these wake elements according to the 
geometry shown in the figure, but in calculating the contribution of this 
element to the lift vorticity, the code would take the trailing edge of this 
last wake element to be at X ~ 00 * 

Flow Inlets 

We have added a feature to the Hess code to allow simulation of flow 
into the aperture of an inlet. (The code is not organized to handle internal 
f 1 ows . ) 

The aperture is represented by quadrilateral panels in the same manner 
as the body surface. To illustrate this. Fig. 6 shows the paneling of the 
orifice in the tip of the intake tube of a cloud water meter, the EWER, 
which is mounted under the wing of a C130 research airplane (ref. 14). In- 
let aperture panel coordinates must be the first in the deck of surface point 
cards (card 12). 

Input card no. 2 contains the number of aperture quadrilaterals (LEAK) 
and also the fraction of the free-stream flow speed that is "leaked" through 
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the apertures (FRAC). This leakage 
is taken to be the same for each 
aperture quadrilateral. If there 
is no flow inlet, the LEAK and FRAC 
fields on card 2 are blank. 

Off-Body Points 

If logical parameter LOFF (card 
3) is true, flow velocities are com- 
puted at user-specified arbitrary 
points (cards 15) and printed near 
the end of the output. The computa- 
tion and processing of these velocities 
have no effect on the calculation and 
storage of data to be used by FLOVEL. 


Figure 6. 


Computer plot of tip and 
orifice of a EWER cloud 
water content meter probe 
(ref. 14). 


Symmetry Planes 

Two symmetry planes are allowed 
by this code. These correspond to 
the first and second symmetry planes 
described above for PGM BPOXC (p. 20) and are applied in the same order. 

That is, the first symmetry plane, or the single symmetry plane if only one 
is specified, is the y = 0 plane. The second symmetry plane is the z = 0 
plane. When both symmetry planes are specified, surface point coordinates 
are reflected across the y = 0 plane first, then across the z = 0 plane and 
finally back across the y = 0 plane. Only the primary surface point data 
should be input; if reflected as well as primary data are input, and both 
are reflected, the flow calculation results will be incorrect. 
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Each symmetry plane is designated as either "plus" or minus". A 
plus symmetry plane is an ordinary reflection plane, on which zero normal 
velocity is defined at all points. Velocity potential is zero at all 
points on a minus symmetry plane; the minus symmetry plane reflects coordi- 
nates normally, but reverses the direction of induced velocities computed 
from the reflected point. An example of a minus symmetry plane is a free 
surface for the condition of infinite Froude number. 

When symmetry is specified for a body with one or more lifting sections 
(airfoil), the code assumes that the y = 0 plane is the airfoil symmetry 
plane, and that the y direction is the airfoil spanwise direction. This 
requires reversal of direction of vorticity induced velocities computed 
from points reflected across the y = 0 plane. (Note that this reversal 
would be independent of that caused by specification of the y = 0 plane 
as a minus symmetry plane; indeed the two reversals would cancel.) Vorticity 
induced velocities computed from points reflected across the z = 0 plane are 
not reversed in direction (unless, of course, the z = 0 plane is designated 
a minus symmetry plane). 

Read the next section before a decision is made to use symmetry. 

Uniform Onset Flow (Free Stream Flow) 

In the context of the unmodified Hess code, specification of a uni- 
form onset (i.e., free stream) flow is equivalent to specification of 
airplane angle of attack. Since extra angles of attack add relatively 
little to the computation burden, the code allows up to ten specifications. 

However, if symmetry is specified (see above) free stream vectors must be 
parallel to symmetry planes. 

For calculation of particle trajectories certain restrictions are im- 
posed. First, only the first specified of the angles of attack is used 
for trajectory calculations. Second, since gravity is involved in the 
trajectory calculations, the z axis is always directed vertically upward. 

The x and y axes are in the horizontal plane with the positive x axis 
directed along, or at least close to the direction of, the uniform onset 
flow vector. The y axis is oriented such as to define a right-handed car- 
tesian system of axes. More detail is given above on pp. 19-20. 
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Uniform onset flows are defined in terms of the direction cosines of 
their vectors with respect to the above coordinate system. Normally, 
with the uniform onset flow vector coincident with the x axis vector, the 
direction cosines are: 1.0, 0.0, 0.0, and airplane attitude angle is 

adjusted by rotation of coordinates about the y axis as specified by input 
data on card no. 5. Since uniform onset flow velocities always are taken 
to have unit magnitude, the direction cosines are equivalent to the free 
stream velocity components. 

Correction for Compressibility Effects 

Since the Hess method computes potential flow it is not capable of 
directly coping with compressibility effects, so that an auxiliary correction 
must be used. The Prandtl -Glauert correction (ref. 15) is used here. 

The correction is applied by user specification of a nonzero free 
stream Mach number (MACH) on card 2. While the correction is automatically 
applied for any free stream Mach number, N^, in the range 0< N M < 1, it 
is not significant for N M < 0.5, so that the MACH field of card 2 may as 
well be left blank unless N M > 0.5. 

The small perturbation approximation is basic to the derivation of 
the correction (ref. 15), with the consequence that only streamwise com- 
pressibility effects are accounted for. Thus, compressibility effects that 
would result from locally high Mach numbers are ignored. 

This correction is included in the nonlifting code described in refer- 
ence 2, but it is not described by Hess and Smith in reference 4, nor is it 
included in Hess' lifting code (refs. 3 and 6). Since the correction is 
not described in connection with these codes elsewhere, details are given 
here in Appendix A. 

The Kutta Condition: Pressure Equality and Flow Tangency Options 

Since the air that leaves the trailing edge of a three-dimensional airfoil 
after passage over its upper surface differs in flow from that which has 
passed under, there is a flow discontinuity extending aft of the trailing edge 
which constitutes the wake vortex sheet. Specification of flow at the trail- 
ing edge, the so-called Kutta condition, is fundamentally important in 
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specifying the overall lifting flow. In a method which uses finite elements, 
such as this one, the discontinuity is not abrupt at the trailing edge, and 
selection and application of an appropriate Kutta condition is less obvious 
than for analytical methods. This situation is discussed in depth in sec. 

6.5 of reference 3, to which the reader is referred for details. 

Two options are available for specification of the Kutta condition: 

1. Flow tangency. This option requires specification of the point 
of application of the Kutta condition for each lifting strip, 
and for each point must be specified the normal unit vector to 
the wake. (These data are not required for extra strips.) 

Results are sensitive to these quantities. 

2. Pressure equality. For this option, surface pressures at the 
centroids of the two on-body elements in each lifting strip 
which include the trailing edge are constrained to be equal. 

This constraint is imposed automatically by the code so that 
no additional information or data are required. Moreover, it 
turns out that this method is insensitive to the distance of 
the points of application from the trailing edge. 

The user is encouraged to select the pressure equality option (ref. 3, 
sec. 8.2). 

Spanwise Variation of Vorticity: Step Function and Piecewise Linear Options 

As illustrated in Figures 1 and 4, a lifting surface is partitioned by 
n-lines into strips of elements which are called lifting strips. Each 
lifting strip consists of all of the elements in a complete circuit around 
the airfoil, including the wake. A dipole distribution, from which vorticity 
is calculated, varies systematically around the lifting strip circuit such 
that all individual element contributions in a strip may be summed. This 
is done such as to yield two vortex strengths, B, per unit strip length for 
each strip: one for each edge of each strip. Spanwise variation of these 

vortex strengths is handled by one of two options: 
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1. Step function. The vortex strength is constant spanwise across 
each lifting strip and changes discontinuously at the n-line 
boundaries between lifting strips. 

2. Piecewise linear. Vortex strength varies linearly spanwise 
across each lifting strip such as to reduce the discontinuities 
at the n-lines to higher order effects. 

The step function option is both simpler and less demanding computationally, 
and, in practice, seems to yield results which are not significantly inferior 
to the piecewise linear option (ref. 3, sec. 8.3). Use of the piecewise 
linear option also requires input of special data, among which are widths of 
all of the lifting strips. This is required since a particular strip does 
not necessarily have a constant width along its (chordwise) length, and if 
not, it may not be obvious how to determine the required nominal width for 
the strip. Thus, the user must enter a width for each strip, which for the 
usual case of constant width is the constant width value. Also when the 
piecewise linear option is used, complications arise at the ends of lifting 
sections, which complications are handled by specification of input parameters 
NLINE1 and NLINEN (card 8) as described next. 


End Conditions for Lifting Sections When the Piecewise Linear Opt ion is 
Exercised : Parameters NLINE1 and NLINEN ~ 

Parameters NLINE1 and NLINEN must be specified by the user for each 
lifting section when the piecewise linear option is selected to ensure a 
reasonable treatment of wing tip and root, and to provide continuity between 
contiguous, independent lifting sections. NLINE1 specifies the edge condi- 
tion at the beginning (i.e., at the first lifting strip) of the section, 
and NLINEN specifies the edge condition at the opposite end of the section. 
For the Jth lifting section we define: 
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NLINEl(J) 


Value 

Description 


1 For the first lifting section (J=l) this provides the 
"normal" definition of a wing tip (ref. 3, pp. 92-93). 
An extra strip with width equal to the first regularly 
defined strip is created. Used to define the wing tip 
when the wing is described by strips (and sections) in 
sequence from tip to root. 

2 Never used. Causes an error condition print and job 
abortion. 

3 Used for J > 1 when the first lifting strip of the Jth 
lifting section is contiguous with the last strip of 
the J - 1th lifting section. Creates an extra strip of 
width equal to the width of the last strip of the J-lth 
section. In combination with NLINEN(J-l) =2, it 
ensures spanwise continuity of vorticity between con- 
tiguous lifting sections. 

4 The code accepts input to define an extra first strip. 
The complete extra strip is defined as though it were 
an ordinary strip, via input of no. 12 cards, and its 
width is defined via WIDXTR(1,J) of card 11. Used to 
define the root, inside of the fuselage, of a wing 
described by strips (and sections) in sequence from 
root to tip. 

5 Same as NLINE1=4, but the input (via WIDXTR(1,J) in 
card 11) strip width is ignored, and the extra strip is 
assigned the width value of the first ordinary strip 

in the section. 
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NLINEN(J) 

Value 


Descri ption 


1 For the last lifting section, this provides the "normal" 
definition of a wing tip (ref. 3, pp. 92-93). An extra 
strip with width equal to that of the last regularly 
defined strip is created. Used to define the wing tip 
when the wing is described by strips (and sections) in 
the sequence from root to tip. 

2 Used when the end of the Jth lifting section is contigu- 
ous with the start of the J+lth lifting section. Creates 
an extra strip of width equal to that of the first strip 
of the J+l section. In combination with NLINE1(J+1) = 3, 
it ensures spanwise continuity of vorticity between con- 
tiguous lifting sections. 

3 Never used. Causes an error condition print and job 
abortion. 

4 The code accepts input to define an extra strip. The 
complete extra strip is defined as though it were an 
ordinary strip, via input of no. 12 cards, and its 
width is defined via WIDXTR(2,J) of card 11. Used to 
define the root, inside of a fuselage, of a wing that 
is described by strips (and sections) in the sequence 
from tip to root. 

5 Same as NLINEN=4, but the input (via WIDXTR(2,J) on 
card 11) strip width is ignored, and the extra strip is 
assigned the width of the last ordinary strip of the 
section. 
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Forces and Moments 


As explained in detail in sec. 7.14 of reference 3, pressure force is 
calculated on each element, and the moment of this element is computed 
about an origin that is specified by the user (card no. 6). Total forces 
and moments, for strips, sections and the whole body, are calculated by 
simply adding the individual vectors. 

Variable Array Dimensioning 

To render the code open ended with respect to its capacity to handle 
numbers of quadrilateral elements, and numbers and types (lifting or non- 
lifting) of sections and strips, we have incorporated the variable dimension 
FORTRAN feature into the DUGLFT code. Thus, in all but the main program, 
DUGLFT , array dimensions are represented by symbolic parameters that are 
discussed below and in the COMMENT statements found in the FORTRAN listing 
of DUGLFT. These parameters are assigned numerical values in DUGLFT (cards 
131 and 132). Also in DUGLFT, all of the arrays are dimensioned numerically, 
rather than symbol ically, with values that are identical to those assigned 
to the symbolic parameters. Therefore, to adjust array dimensions to fit 
requirements of a particular problem, it is necessary to change them only 
in program DUGLFT (and in COMMON/ BASDAT/of SR SETFLO and SR FLOVEL if flow 
velocities are to be calculated using the data output from DUGLFT). 

Minimum values for the variable dimension parameters are as follows. 

They never include quantities that are generated by symmetry. 

LFSX Number of lifting sections 

NL2X NSLX + 2 

NOBX Total number of lifting strips, not counting extra strips 

NONX Number of on-body elements in the flow (not counting 

ignored, wake and extra strip elements) plus Kutta 
points defined by input (flow tangency option only) 
plus off body points (cards 15) 
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NSEX Total Number of sections (lifting plus nonlifting) 

NSLX Maximum number of lifting strips in any lifting section 
(including extra strips if input) 

NSTX Total number of strips (i.e., n-lines) (lifting plus 
nonlifting including extra strips if input) 

NWAX Similar to 3 x (number of on-body quadrilateral elements in 

the flow (see NONX) plus the number of onset flows*) 

and 

Cube of the number of lifting strips (not counting extra 
strips) (i.e. NOBX 3 ) if the pressure equality Kutta con- 
dition option is selected. 

N2BX 2 x NOBX 

The DUGLFT codes are the only ones in the package described here that 
use the variable dimension FORTRAN feature. The trajectory codes and pro- 
gram FLOPNT are not explicitly affected by the DUGLFT dimensioning, but 
since they all call subroutines SETFLO and FLOVEL, they are implicitly 
affected. 

Subroutine SETFLO recovers data from external unit 18 for use by 
FLOVEL that was previously stored there by the DUGLFT codes. Therefore, 
the dimensions of the arrays in COMMON/ BASDAT/ in subroutines SETFLO and 
FLOVEL must be consistent with those defined in DUGLFT; and if the DUGLFT 
array dimensions are changed, the user must be careful to also change 
corresponding dimensions in COMMON/BASDAT/, or at least to be sure that 
the COMMON/BASDAT/ array dimensions are sufficient to contain the data on 
unit 18. 


*The number of onset flows is the number of lifting strips (not counting 
extra strips) plus the number of uniform onset flows. 
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A warning is in order regarding expansion of the arrays to handle 
very large problems. The user should remember that flow velocity computa- 
tion time increases directly with the number of quadrilateral elements 
used. Moreover, a large number of calls to subroutine FLOVEL is required 
for each trajectory calculation, and FLOVEL may use all of its 
subroutines or each call, some of them repeatedly. Therefore, it is not 
practical to use segmentation or overlays with FLOVEL so that all of the 
data and all of the code must be carried in rapid access memory all of the 
time. 

Obviously, large problems can be costly in terms of both computing 
time and storage, and accordingly, detail in a surface description that is 
not required for a particular problem should be ignored or be replaced by 
paneling that is as coarsely resolved as is practicable. This may require 
that more than one digital description of a body surface be prepared. An 
example is shown in Figure 7, which shows digital descriptions of the De 
Havilland Twin Otter airplane at two levels of resolution. They have been 
prepared in mutually compatible sectional form such that various sections 
can be interchanged. For example, if attention is focused on the forward 
fuselage, the detailed nose and cabin sections from Figure 7a can be combined 
with the coarsely resolved cylindrical fuselage of Figure 7b, and the wings 
and tail assemblies can be ignored completely. 

Printed Output 

Printed output consists of two main parts, plus a summary of input 
control data, various error condition messages, and optional outputs of 
data that are used for debugging. An example printout, which is for the 
test problem, is included in the Example Problem section below. 

1. The first printout is a summary of input control data, and is 
self-explanatory. 

2. Next, which is the first main part, is a printout (from NOLIFT 
and LIFT) of element data which is almost the same as the printout 
for PGM PBOXC discussed above on pp. 30-31. The major differences 
reflect the fact that DUGLFT code uses element centroids as 
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a. 


Finely resolved version 



Coarsely resolved version 


Figure 7. Digital description of the DeHavilland Twin Otter 
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control points rather than null points, which eliminates need 
for the print of the integers 1 or 2 to label abnormal null 
points, and the control point coordinates are indicated as XO, 

YO, ZO rather than NPX, NPY, NPZ. Also elements are labeled as 
nonlifting, lifting or wake. For lifting sections, D is the 
denominator of equation (7.2.3) of reference 3. 

A short table follows (from INPUTL) titled TABLE OF INPUT INFORMA- 
TION, which summaries the data in terms of section type, number 
of elements per section, number of strips, etc. 

3. In the course of computing velocities induced by each element 
on all others, additional summary information is printed (from 
VIJMX) for each section. For lifting strips, this includes 
information on ignored elements which does not appear elsewhere. 

4. If the piecewise linear option for determination of spanwise 
variation of vortex strength is used, strip widths, W^, and 
parameters D^, E^, F^ (ref. 3, sec. 7.11) are printed for each 
strip for each section, along with a summary of edge conditions 
(NLINE1 and NLINEN, card 8 and pp. 40-42). 

5. A short statement is printed (from C0LS0L) regarding the dimen- 
sions of the matrix that is solved to determine element source 
strengths, and the number of right-hand-sides (i.e., number of 
uniform onset flows plus number of lifting strips) for which the 
solutions are obtained. 

6. The second main printout (from PRINTL) contains the final results 
of the calculations. A printout for each on-body element is 
labeled as follows: 
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XO, YO, ZO 
VX, VY, VZ 

VT 

VTSQ 

CP 

DCX, DCY , DCZ 

NX, NY, NZ 

SIG 

VN 

AREA 


Control point coordinates 

Flow velocity components at the 
control point 

Velocity magnitude 

Square of velocity magnitude 

Pressure coefficient = 1.0 - VTSQ 

Direction cosines of the velocity 
components 

Components of the unit normal to 
the plane of the element 

Source density 

Velocity component in the direction 
of the unit normal 

Area of the element. 


Printouts for off-body and Kutta points are similarly labeled. 

Also printed are vector components for pressure force and 

moment for each strip, each section and for the entire body, 

( k ) 

as well as a table of vortex strengths per unit length, B , 
for each lifting strip. 


7. Error Messages (ref. 6) 

(a) Message: MISMATCH OF ELEMENTS IN A LIFTING STRIP IS DETECTED. 

ELEMENTS FORMED = xxx, ELEMENTS INPUT = xxx, COMPU- 
TATION TERMINATED. (SR INPUTL) 

Cause of error: Inconsistent input data. The program sums 

the number of on-body elements plus the wake elements 
specified on card 8. This sum does not match with the 
elements formed from the input coordinates. 

Action: Check the lifting body information card (card 8) and 

the quadrilateral corner point coordinate cards 
(cards 12). The number of points on an n-line should 
equal the number of elements plus 1. 

For example. If in a lifting section each lifting 
strip consists of 10 on-body elements and 1 wake 
element, the total number of elements is 11, and 
there should be 12 points on each n-line input via 
cards no. 12. 

(b) Message: ERROR IN IGNORED ELEMENT COUNT xxx, SHOULD BE xxx. 

(SR LIFT) 

Cause of error: Erroneous specification of the ignored element 

information. 

Action: Check card 10 to make sure the ignored element informa- 

tion is properly specified. 

(c) Message: LABEL ERROR IN NONLIFTING VFORM. (SR VFMNLF ) 

LABEL ERROR IN LIFTING VFORM. (SR VFMLFT) 
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Cause of error: Geometric data for each element strip, 

preceded by a lifting or nonlifting label, are 
stored on unit 4. The error occurs when a label- 
ing mixup is detected during input of the data from 
unit 4 for calculation of velocities. That is, 
data for a strip labeled lifting are encountered 
during computation for a nonlifting section, or 
vice versa. 

Action: Check that the number of lifting strips specified 

on card no. 8 for each lifting section corresponds 
with the cards no. 12 input. 

(d) The following messages pertain to errors in specification of 
variable dimensions (SR CKARRY): 

ELEMENT CAPACITY, NONX = xxxxxx IS LESS THAN TOTAL ELEMENTS, 
NON = xxxxxx 

STRIP CAPACITY, NSTX = xxxx IS LESS THAN TOTAL STRIPS, NSTRP 
= xxxxx 

SECTION CAPACITY, NSEX = xxxx IS LESS THAN TOTAL SECTIONS, 
ISECT = xxxx 

LIFTING SECTION CAPACITY, LFSX = xxx IS LESS THAN NUMBER OF 
LIFTING SECTIONS, LIFSEC = xxxx 

LIFTING STRIP CAPACITY, NOBX = xxx IS LESS THAN TOTAL LIFTING 
STRIPS, LSTRP = xxxx 

N2BX = xxx IS NOT .GE. TWICE NOBX AS REQUIRED, NOBX = xxx 

NSLX = xxxx IS LESS THAN THE MAX. NO. OF STRIPS IN A LIFTING 
SECTION, WHICH IS xxxx 

CAPACITY OF ARRAY WKAREA, NWAX = xxxxxx, USED BY COLSOL TO 
DETERMINE SOURCE STRENGTHS, IS INSUFFICIENT. IT MUST BE 
•GE. xxxxxxx 

NWAX = xxx IS NOT .GE. NO. OF LIFTING STRIPS = xxx CUBED, AS 
REQUIRED FOR THE PRESSURE EQUALITY KUTTA OPTION 



Cause of error: Array dimensions are inadequate to accomodate 

the input data. 

Action: Check array dimensions and variable array parameter 

(see pp. 43-45 and PGM DUGLFT FORTRAN listing) 
against the storage demands of the element data 
input via cards no. 12. Also check input parame- 
ter LIFSEC, NSORCE , NWAKE , NSTRIP, and IXFLAG. 

(e) Messages: xxx ANGLES OF ATTACK HAVE BEEN SPECIFIED, ONLY 

ONE IS ALLOWED SINCE COMPRESSION EFFECTS ARE CON- 
SIDERED. 

ANGLE OF ATTACK, +x.xxxxxE+xx, +x.xxxxxE+xx, +x.xxxxx 
E+xx IS INAPPROPRIATE FOR A CASE WITH COMPRESSION 
CORRECTION. (SR CKARRY) 

Cause of error: Only one uniform onset flow (i.e., free 

stream) is allowed if the compressibility correction 
is applied (MACH > 0.0 on card no. 2). Moreover, 
the direction cosines (ALPHAX, ALPHAY, ALPHAZ) of 
this onset flow must be (1.0, 0.0, 0.0)(card no. 4). 

Action: Set IATACK = 1 on card 2, and/or specify direction 

cosines on card 4 as stated above. 

(f) Message: THE NUMBER OF KUTTA POINTS SPECIFIED = xxx IS IN- 

CORRECT AND SHOULD BE xxx (SR CKARRY) 

Cause of error: The flow tangency Kutta option has been 

specified, and the number of Kutta points speci- 
fied by input (cards no. 9, 13 and 14) does not 
equal the number of lifting strips. 
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Action: Check parameter KUTTA on card 9, and the number 

of KUTTA data points on cards 13 and 14, against 
the number of lifting strips input via cards no. 

12. (Do not count extra strips.) 

(g) Message: ERROR IN VFORM. THE ELEMENTS FORMED DO NOT CORRES- 

POND TO THE NO. OF BODY ELEMENTS. (SRS VFMNLF AND 
VFMLFT) 

Cause of error: Element tally recorded by SR INPUTL does 

not match with tally recorded from input of data 
from unit 4 during velocity calculation. 

Action: Check lifting strip specification data on card 8 

for consistency with cards no. 12 input data. 

(h) Message: AFTER xxx ITERATIONS, DELTA B STILL DID NOT CON- 

VERGE TO THE GIVEN CRITERION / LARGEST DELTA B = 

+ x. xxxxxxE+xx / PROGRAM PROCEEDS WITH THE MOST 
CURRENT VORTEX STRENGTH. (SR P KUTTA) 

Cause of error: Nonconvergence of vortex strengths, B, calcu- 

lation via the pressure equality Kutta condition 
method (ref. 3, sec. 7.13.2). 

Action: Check the cards no. 12 input data. 

(i) Message: THIS CODE SHOULD BE APPLIED TO FIRST STRIP. 

or 

THIS CODE SHOULD BE APPLIED TO LAST STRIP. (SR DKEKFK 
or PSONST) 

Cause of error: Improper specification of NLINE1 or NLINEN 

for piecewise linear option. Specifically, either 



NLINE1 = 2 or NLINEN = 3 is specified, both of 
which are forbidden, (see pp. 40-42) 

Action: Check card 8 specifications. 

(j) Message: xxx ON-BODY POINTS MISSED. EXECUTION TERMINATED. 

(SR PRINTL) 

Cause of error: The number of on-body source elements tallied 

during final printout does not agree with the count 
tallied during input. 

Action: Check input data. 

(k) Message: xxx KUTTA POINTS MISSED. EXECUTION TERMINATED. 

(SR PRINTL) 

Cause of error: The number of Kutta points tallied during 

final print out does not agree with the number 
specified by parameter KUTTA on card 9. 

Action: Check the number of Kutta points input via cards 

13 and 14 against parameter KUTTA. 

(l) Message: xxx OFF-BODY POINTS MISSED. EXECUTION TERMINATED. 

(SR PRINTL) 


Cause of error: The number of off-body points tallied during 

final printout does not agree with the number tallied 
during input. 

Action: Check input data. 
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8. Optional Printouts for Use In Debugging 

(a) Geometrical data for each element. (IOUT = .TRUE., card 3) 

(SR INPUTL) For each nonlifting element is printed the element 
sequence number and twenty-nine geometric quantities (ref. 4, 
sec. 9.51), and for each lifting element is printed the element 
sequence number and forty-five geometric quantities (ref. 3, 
sec. 7.2). 

(b) Source induced velocity matrix, V^. (MPR = 1, card 2) (SR 
PNTVIJ ) 

COLUMN Matrix column number (j) 

CNTRL PT Control point number (i) 

VXS, VYS, VZS Velocity components 

If LIFSEC.GT.O (card 2), dipole induced velocity matrices, 
vtf>. vjj^, also are printed. 

STRIP Lifting strip number 

CONTRL PT Control point number 

VXF , VYF, VZF 

First and second velocity components 

VXS, VYS, VZS 

(c) Onset flow matrices, v( k ^ and v|“^. (MPR = 3) (SR UNIFLO) 

ONSET FLOW NO. 

CONTROL POINTS Control point numbers 

X - FLOW 

Y - FLOW Onset flow velocity components 

z - FLOW 
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(d) Dot product matrices, A^, n|^ and N^. (MPR >_ 2, card 2) 
(SR AIJMX and N I KMX ) (See ref. 3, eq. (7.12.5)) 

COLUMN Matrix column number (j) 

AIJ Elements of A.. 

' J 

FLOW NO. Onset flow number (k) 

RIJ Right side of eq. (7.12.5), ref. 3 

(e) Source density matrix (MPR >_ 2, card 2) (PGM SIGMAL) SOLU- 

TION OBTAINED AFTER COLSOL FLOW NO. Onset flow number. 
Element source densities, and are printed eight 

to a line. 


Unit 18 Output 

The following data are stored on peripheral storage unit 18 (in binary 
format) for use later by SR FLOVEL in calculating flow velocities at arbitrary 
space points. Actual record structures are most easily determined by examin- 
ing the SR SETFLO FORTRAN listing. 


SR ST0R18: 

CASE 

ISECT 

LIFSEC 

ALPHAX(l) 

ALPHAV(l) 

ALPHAZ(l) 

SYM1 

SYM2 

NSYM 

NSTRP 


Body identifier (input card 2) 

Number of sections (lifting plus nonlifting) 
Number of lifting sections (input card 2) 

Uniform onset flow direction cosines (card 4) 

Floating point equivalents of input parameters 
NSYM1 and NSYM2 (card 2) 

Total number of symmetry planes 

Total number of strips, including extra lifting 
strips if input. 
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For each nonlifting element, the twenty-nine geometric quantities 
written on unit 4 by SR NOLIFT. 

For each lifting element, the forty-five geometric quantities written 
on unit 4 by SR LIFT. 
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SR DKEKFK: 


Only if the piecewise linear method is used for calculation of 
spanwise variation of vorticity. For each of K strips in J = 
LIFSEC lifting sections, 

K,(D(I,J),E(I,J),F(I,J),I=1,K) 

where D, E, F are D^, E^, F^ of eq. (7.11.5) of reference 3. 


SR SUMSIG: 

KFLOW Number of lifting strips 

KONTRL Number of on-body source elements (not including 

ignored, wake and extra strip elements) 

C0MSIG( KONTRL) Combined source densities (ref. 3, eq. (7.13.1)) 
B(KFLOW) Vortex strength per unit length 

Segmentation Structure 

The DUGLFT code may consume a large amount of storage, and it is best 
to use overlays or segmentation to reduce this. Figure 8 shows a segmenta- 
tion tree structure that can be used, and Table 3 shows the corresponding 

tree directions required by the Control Data Corporation CYBER loader. The 
exact tree structure can be derived from these displays. 

TABLE 3 

CDC CYBER LOADER TREE DIRECTIVES REQUIRED 
FOR THE SEGMENTATION STRUCTURE OF FIGURE 8 

Directive . . 


Number 

Label 

Verb 


Specification 

1 

INPUT 

TREE 

INPUTL - 

(NOLIFT, LIFT) 

2 

VFLIFT 

TREE 

VFMLFT - 

(NEAR,WNEAR) 

3 

VMATX2 

TREE 

VIJMX - 

(VFMNLF, VFLIFT) 

4 

PIECE 

TREE 

PISWIS - 

(DKEKFK, PSONST) 

5 

VMATX1 

TREE 

VMATRX - 

(VMATX2, PNTVIJ.STEPFN, PIECE, UNIFLO) 

6 

BVORXX 

TREE 

BVORTX - 

(PKUTTA,FKUTTA) 

7 

DUGROT 

TREE 

DUGLFT - 
SIGMAL 

CONTRL - ( INPUT, CKARRY, VMATX1 , 
, BVORXX, SUMSIG, VELOCY) 
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CONTRL 
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Figure 8. DUGLFT segmentation tree structure. 





























Peripheral Storage 


In addition to the system input and print units (5 and 6), and unit 18 
which is used for output of data required by subroutines SETFLO and FLOVEL 

(see pp. 55-57), DUGLFT uses eleven units for scratch storage. All data stored 
on these units are in binary format. In the following, use of each unit 
is considered only in terms of the maximum number of data words (numbers) and 
record lengths that would be stored on it. The following variables are de- 
fined to aid in this: 


KONTRL 


KUTTA 


NOFF 

NON 

IATACK 

KFLOW 

NFLOW 


Number of quadrilateral elements, not including 
those: generated by symmetry, ignored, in the 

wake and in extra strips. 

Points defined by input cards no. 13 and 14 at 
which the Kutta condition is to be applied. 

(KUTTA > 0 only if the flow tangency option 
is exercised.) 

Number of off-body points at which velocity is 
to be calculated as defined by input cards no. 15. 

KONTRL + KUTTA + NOFF 

Number of uniform onset flows 

Number of lifting strips, not counting extra 
strips, nor those generated by symmetry 

KFLOW + IATACK 


Unit 3: NFLOW records each consisting of 3 x NON numbers 

Unit 4: There is a record of 29 numbers for each nonlifting quadrilateral 
element 

plus 

There is a record of 45 numbers for each lifting quadrilateral 
element (including ignored, wake and extra strip elements) 

plus 

A one-word record for each section of elements 
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Unit 8: The larger of 

Two records each of length 3 x NON numbers 

or 

NFLOW records each of length 6 x KFLOW numbers 

or 

KONTRL records each of length KUTTA numbers 
Units 9 and 10: KONTRL records of maximum length KONTRL + 1 numbers 

Units 11 and 12: % KONTRL records each of length 3 x NON numbers 

Unit 13: The larger of 

2 x KONTRL records each of length 3 x NON numbers 

or 

KONTRL records of maximum length KONTRL + NFLOW + 1 numbers 

Unit 14: The larger of 

KFLOW records each of length 3 x NON numbers 

or 

KONTRL records of maximum length KONTRL + 1 numbers 

Unit 15: The larger of 

KFLOW records each of length 3 x NON numbers 

or 

KONTRL records of maximum length KONTRL + NFLOW + 1 numbers 
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DUGLFT Card Input 


Card 

N°» Variables and Format 

1 (TITLE(I) ,1=1 ,18) , 

( 1 8A4 ) 

2 CASE , LIFSEC, IATACK, 
NSYMI ,NSYM2,MPR, 
LEAK, FRAC, MACH 
(A4,6I4,2X,2F10.0) 


3 IPROS.LOFF, MOMENT, LIST, 

IOUT ( 5L1 ) 


• — — Description 

Run identification 


Run control data: 


CASE 

(Col. 1-4) 

LIFSEC 
(Col. 5-8) 


IATACK 
(Col. 9-12) 


Four character body idenfitication. 
Total number of lifting sections. 


c r;:r m ~ u.e., uniform free 

stream flows) to be specified via cards no. 4. 
aximum value is 10. If the compression correc- 
tion is to be applied (MACH > 0.0), it is 
necessary that IATACK =1 


NSYM1 

(Col. 13-16) 
NSYM2 

(Col. 17-20) 


m , • Wee values 0, 1 or -1 is entered. 

" specifies the 1st symmetry plane, and 
NS YM2 specifies the second symmetry plane 
according to: 


0 nonexistent 


+1 a plus (ordinary) symmetry plane 
-1 a minus (anti) symmetry plane 
See pp. 36-37. 


”.21-24) S n e *pp! , |4-K d ) , ° r Pr ° Sra " m ' y 

0 No debug print. This is the normal value 
tor this parameter. 


1 Print the source induced velocity matrix, 

V ij, and > if LIFSEC > 0, print the dipole 

induced velocity matrices, v(P and 

1 * 1 j 

>2 Print the dot product matrices A.., 

00 (oo) 1J 

and N: 00 ', and the element source den- 
sities aj k and a(°°\ 

3 Print the onset flow matrices, v( k ^ and vW 

LEAK Number of inlet quadrilateral elements. These 

v Col. 25-28) must be the first elements in the digital descrip- 
tion set (cards no. 12). See pp. 35-36. 

FRAC Fraction of the unit free stream flow that passes 

(Col. 31-40) through each of the LEAK inlet elements. If 
LEAK = 0, leave this field blank. 


MACH Mach number of the free stream flow. (Note 

(Col. 41-50) that this is a floating point number.) If 

% 1 0*5, leave this field blank. See p. 38. 

Logical control flags: 

IPROS If true, the card 12, 13 and 15 coordinates are 

(Col.l) to be translated, scaled, and rotated about 

the y axis before processing, and card no. 

5 is to be input. If false, no translation, 
scaling or rotation is done, and card no. 5 
is not input. 
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DUGLFT Card Input, con't. 


Card 

No. Variables and Format 


4 ( ALPHAX ( I ) ,ALPHAY ( I ) , 

ALPHAZ(I),I=1,IATACK) 
(3E10.0) 


5 ANGLE, XSCALE,YSCALE, 
ZSCALE ,XTRANS ,YTRANS , 
ZTRANS , (7F1 0. 0) 

6 ORIGNX ,ORIGNY ,0RIGNZ, 
(3E10.0) 


Description 


LQFF If true, velocities at off-body points are to be 

(Col. 2) calculated. The off-body points are specified 

by the user via input of the no. 15 cards. If 
false, off-body velocities are not calculated 
and there is no input of the no. 15 cards. 

MOMENT If true, the moment origin is specified by input 

(Col. 3) of card no. 6. If false, card 6 is not input 

and moments are computed about point (0,0,0). 

LIST If false, specifies complete execution. If true, 

(Col 4) execution is terminated after the first main 

part of the printed output (i.e., after comple- 
tion of printed output item 2 as described on 
pp. 45 and 47). 

j OUT If true, the 29 geometric quantities for each 

(Col. 5) nonlifting element and 45 geometric quantities 

for each lifting element are printed. The normal 
value for this parameter is false. 

Direction cosines of uniform onset (i.e., free stream) flow 
vectors. IATACK is the number of uniform onset flows speci- 
fied in card no. 2. One set of direction cosines per card. 

If the compression correction is applied (MACH > 0.0), only 
one uniform onset flow vector can be specified. If more than 
one vector is specified, only the first is passed along via 
unit 18 for use by SR SETFL0 and FL0VEL. The direction cosines 
are with reference to the airplane coordinate system (after 
rotation by angle ANGLE (card 5)). These vectors are equiva- 
lent to unit free stream velocities. Ordinarily, free stream 
unit velocity components are (1 .0,0. 0,0.0) . See pp. 37-38. 

Input only if IPR0S = .TRUE, on card 3, Same as card no. 4 
of program PB0XC. 


Coordinates of the moment origin. This card is input only 
if MOMENT = .TRUE, on card 3. 


Cards 7 through 11 are input only if LIFSEC > 0 on card 2. 


7 LKUTT ,LASWAK,PESWIS , 

IGW(5L1 ) 


Logical control flags for lifting section data: 

LKUTT If true, the flow tangency method for applica- 

(Col.l) tion of the Kutta condition is selected. This 

means that one point in or near the wake of 
each lifting strip (not counting extra strips) 
must be specified via input of the no. 13 and 
14 cards, and the number of these points must 
be specified via input of card no. 9. If false, 
the pressure equality method is selected, and 
cards no. 9, 13 and 14 are not input. 


LASWAK 
(Col. 2) 


If true, the trailing edge of the last wake 
element is automatically extended by the code 
t 0 x = oo. This is the semi-infinite last wake 
element option. (See p. 35.) 
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DUGLFT Card Input, Con't. 


Card 

No. Variables and Format 


8 (NSORCE(J) ,NWAKE(J) , 

NSTRIP(J) ,NLINE1 (J) , 
NLINEN(J) , IXFLAG(J) , 
J=1 ,LIFSEC) , (614) 


9 KUTTA( 14) 


PESWIS 
(Col. 3) 


I GW 

(Col. 4) 


NSORCE(J) 
(Col. 1-4) 

NWAKE(J) 
(Col. 5-8) 


NSTRIP(J) 
(Col. 9-12) 


NLINEl(J) 
(Col. 13-16) 


NLINEN(J) 
(Col. 17-20) 

IXFLAG(J) 
(Col .21-24) 


Des cription 

If true, the piecewise linear method for calcula- 
ting spanwise variation of lift vorticity is 
selected, and lifting strip widths must be input 
via cards no. 11. If false, the step function 
option is selected, and cards no. 11 are not 
input. See pp. 39-40. 

If true, there are ignored lifting elements which 
must be defined via input of the no. 10 cards. 

If false, there are no ignored elements, and 
cards no. 10 are not input. See p. 27. 

Number of on-body elements (including ignored) 
in each lifting strip of the Jth lifting section. 

Number of wake elements in each lifting strip of 
the Jth lifting section, including a semi-infinite 
final element if this option is selected. 

(See pp. 26 and 35. ) 

Number of lifting strips in the Jth lifting 
section. Include extra strips only if they 
are defined via input of cards no. 12. 

If the piecewise linear option is selected; 

(PESWIS= .TRUE, on card 7), NLINEl(J) specifies 
the edge condition of the first strip on the 
Jth lifting section (pp. 41-42). If the step 
function option is specified, ignore this field. 

Same as NLINEl(J) but for the last strip of the 
Jth lifting section. 

IXFLAG(J)=0 means that no extra strips are de- 
fined via input for the Jth lifting section. For 
extra strips defined via input (i.e., via cards 
no. 12) , we have: 

IXFLAG(J)=1 means the first strip is an extra 
strip. If the piecewise linear option is selected, 
(PESWIS= .TRUE. on card 7), this also requires 
NLINEl(J) = 4 or 5. 

IXFLAG(J)=3 means the last strip is an extra strip. 
If the piecewise linear option is selected, this 
also requires NLINEN(J) = 4 or 5. 

IXFLAG(J)=2 means that both first and last strips 
are extra strips. If the piecewise linear option 
is specified, this requires that both NLINEI(J) 
and NLINEN(J) = 4 or 5. 


LIFSEC is specified on card 2. A separate card is required 
for each lifting section, and the cards are input in the same 
order as is input of the quadrilateral data via cards no. 12. 

Input only if LKUTT= .TRUE . on card 7. Number of points at 
which the flow tangency method for application of the Kutta 
condition is to be applied. It is required that KUTTA equal 
the total number of lifting strips, not counting extra strips. 
KUTTA is used to read the point coordinates, and the unit 
vectors normal to the wake or airfoil surface at these points, 
via cards no. 13 and 14. 
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DUGLFT Card Input 


Card 

No. Variables and Format 

10 ( (IG1 ( I ,J) ,IGN( I ,J) , 

1=1 , NSTRI P ( J ) ) ,0=1 , 
LIFSEC) ,(1214) 


11 (WIDXTR(I,J),(WIDTH(I,J), 

1=2 , NSTRI P ( J ) - K) , 
WIDXTR(2,J),J=1, LIFSEC), 
(7E10.0) 

K = 0 if IXFLAG(J)=0 
K = 1 if I XFLAG ( J ) = 1 or 3 
K = 2 if I X FLAG ( J ) =2 


12 X,Y ,Z,STAT,LAB ,XX,YY ,ZZ, 

ST ATT, LAB L, 

(3E10.0.2I2/3E10. 0,212) 


Description 


Input only if IGW=.TRUE. on card 7. I=lifting strip index; 
J=lifting section index. If on the Ith strip of the Jth 
lifting section there is a substrip of ignored elements, 
the substrip is defined by specifying its beginning and 
ending element indices (see footnote below) via: 

161(1, J) = index of the first ignored element on the lifting 
strip 

IGN( I , J ) = index of the last ignored element on the lifting 
strip 

If there are no ignored elements on a strip, leave both fields 
blank; but IG1 and IGN must be specified for every lifting 
strip if I GW = .TRUE, on card 7. Six strips per card. Each 
lifting section begins a new card. 

LIFSEC is specified on card 2, and NSTRIP(J) on card 8. 


Input only if PESWIS= .TRUE . on card 7. These quantities are 
the widths of each lifting strip for use in calculating 
spanwise variation of vorticity via the piecewise linear 
method. 

WIDXTR(l.J) specifies the width of the first extra strip of 
the Jth lifting section. If NLINE1(J)^4, leave 
this field blank. 

WIDTH(I.J) specifies the width of the Ith lifting strip of 
the Jth lifting section. 

WIDXTR(2,J) specifies the width of the last extra strip of 
the Jth lifting section. If NLINEN(J)/4, leave 
this field blank. 

LIFSEC is specified on card 2, and NSTRIP(J), NLINEl(J), 
NLINEN(J) and IXFLAG(J) are specified on card 8. 


On-body and wake quadrilateral element corner point coordinates 
are specified by these cards for both lifting and nonlifting 
sections, one point per card. The body and wake surface panels 
are constructed from these data. 


X,Y,Z 
XX ,YY ,ZZ 
(Col. 1-30) 

STAT 
STATT 
(Col. 32) 


Quadrilateral corner point coordinates. 

Status parameter: Allowed values are 0,1, 2, 3: 

0 This point is on the same n-line as the last point. 

1 This point starts a new n-line. 

2 This point starts a new section. 

3 This is the last point in the card 12 input. 


LAB Specifies a lifting or nonlifting section: 

(Col. 34) 0 n0nlifting 

1 lifting 

This field is relevant only when STAT or STATT 
= 2; that is, only on the first card of a new 
section. 

Note: There must be an even number of no. 12 cards; add a blank 

card to the end of the card 12 deck if necessary. 
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DUGLFT Card Input 


Card 

No. 


13 


14 


15 


Footnote: 


Variables and Format 


(XC(I),YC(I),ZC(I), 
1=1 , KUTTA) (3E1 0.0) 


(XN( I ) ,YN( I ) ,ZN( I ) , 
1=1 ,KUTTA) (3E10.0) 


X0F,Y0F,Z0F,STAT,X0FF, 
Y0FF,Z0FF,STATT, 
(3E10,0,I2/3E1 0.0, 12) 


Description 

Input only if LKUTT=.TRUE. on card 7. KUTTA is specified 
via the card 9 input. Coordinates of points (one point per 
card) at which the Kutta condition is to be applied via use 
of the flow tangency method. If IPR0S=.TRUE. (card 3), the 
code automatically translates, scales and rotates these 
coordinates according to the card 5 input data. 

Input only if LKUTT = .TRUE, on card 7. KUTTA is specified 
via the card 9 input. Components of the unit vectors (one 
vector per card) normal to the wake or airfoil surface at 
the points specified by the no. 13 cards at which the Kutta 
condition is to be applied via use of the flow tangency 
method. The order of input must be consistent with that of 
the no. 13 cards. If IPROS = .TRUE, (card 3), a transforma- 
tion is automatically applied by the code to adjust for 
rotation of coordinates by angle ANGLE (card 5). 


X0F,Y0F ,Z0F 
X0FF,Y0FF ,Z0FF 
(Col. 1-30) 


Coordinates of off-body points at which flow 
velocities are to be calculated, one point 
per card. 


STAT 
STATT 
(Col. 32) 


Status parameter. A value of 3 signifies 
the end of the off-body points. Otherwise, 
leave this field blank. 


If IPR0S= .TRUE . (card 3), the code automatically translates, 
scales and rotates these coordinates according to the card 5 
input data. 


Element indices on a lifting strip are simply the sequence numbers of the elements, beginning 
with one at the trailing edge and proceeding along the n-lines in the order with which the 
elements are defined by input. For example, if the 3rd through 7th elements of a strip are 
ignored then IG1 = 3 and IGN = 7. 
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SUBROUTINE FLOVEL 


General Discussion 

Given the coordinates of a point in space (X, Y, Z), the current time 
step interval, H, used in the integration of the particle equations of 
motion, and the particle coordinates at the previous time step, PPR(3), 

SR FLOVEL returns the flow velocity components (VX, V Y, VZ) at point 
(X, Y, Z), and an indicator, INBODY, of whether the body surface has 
been penetrated. INBODY = 0 if the point is exterior to the body, but 
INBODY = 1 if penetration has occurred during the current time step. 

The discussion to follow assumes that a Hess flow field (ref. 

3) is being considered. However, if the user wishes to compute 

flow by use of some other method, for example, flow about a two-dimensional 

airfoil, he may replace FLOVEL by a subroutine of his own design.* 

Implicit in the coding is the assumption that flow velocities are being 
calculated in the course of a particle trajectory calculation. When this 
is not the case, for example when called by FLOPNT, the time step, H, can 
be set to zero in the calling program. This avoids body penetration print- 
outs that may be meaningless or misleading. 

SR FLOVEL is based mainly on the Hess subroutines VFMNLF, VFMLFT, 
NEAR, WNEAR and PSWISE, with modifications and additions required for 
application of combined source densities, vortex strengths, body pene- 
tration tests and special calculations required for points very close to 
element surfaces, element edges or extensions of element edges. The com- 
bined source densities and vortex strengths are used straightforwardly 
as indicated by eqs. (7.13.1) through (7.13.3) of reference 3. 


*In the trajectory codes a call of SR SETFLO precedes the first call of 
FLOVEL. SETFLO reads the data stored on unit 18 by PGM DUGLFT (see p. 

55), which data are required by FLOVEL for calculation of a Hess 

flow velocity, and puts these data into COMMON storage. If a user-designed 

version of FLOVEL is used, SETFLO must be replaced by a dummy subroutine. 
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Body Penetration Tests 


FLOVEL calculates and sums velocities induced at the calculation point 
by each element. Three modes of calculation are used (ref. 3, sec. 7.4): 
(1) for large distances between the calculation point and the element cen- 
troid, point singularity methods are used, (2) for intermediate distances, 
multipole methods are used, and (3) for short distances exact calculations 
are used. For each on-body element for which exact calculations are 

required, a sequence of three sets of body penetration tests are applied. 
The sequence is arranged to save computation time by application first of 
simple, nonexhaustive tests that can determine that penetration cannot 
have occurred; as soon as the possibility of penetration is eliminated, 
further testing is bypassed. In the descriptions of the tests below, we 
assume that a particle trajectory is being computed. 

Sequence 1. If any of the following tests is satisfied, penetration 
is taken not to have occurred: 

1. The vector of separation between the centroid of the element 
and the particle position is projected onto the unit normal 
vector to the plane of the element, and the sign of the pro- 
jection is checked to see if the point is on the exterior 
side of the element. 

2. The square of the distance of the particle position to the 
element centroid is greater than one-forth the square of the 
maximum element diagonal plus the square of the time step. 

3. The magnitude of the projection calculated in test 1 is 
greater than the time step, H. 

For tests 2 and 3 we assume that the maximum particle speed is of order 
unity, so that the maximum distance a particle can travel in one time step 
is of order H. 
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Sequence 2. The vector of separation between the element centroid and the 
particle position at the previous time step is projected onto the unit 
normal vector to the plane of the element. If the sign of this projection 
is the same as that determined in test 1 of the first sequence for the 
current position, or if both previous and current positions are exactly 
in the plane of the element, penetration is taken not to have occurred. 


Sequence 3. The line connecting the current particle position with its 
position at the previous time step is determined, and the point of inter- 
section of this line with the element plane is calculated. If this 
intersection point lies inside the element or on its boundaries, penetra- 
tion is taken to have occurred. Details of this sequence of tests are 
given in Appendix B. 

When penetration is detected, the coordinates of the point on the body 
at which penetration has occurred are printed along with the following 
information: 

REFLECTION LOOP INDEX (12) This tells whether the primary 

element (12 = 1) or a reflected 
element (12 > 1) is involved. 

(See pp. 20 and 36. ) 


ZNP (nonlifting element) 
Z (lifting element) 


Projection of the centroid-to- 
particle position separation vec- 
tor onto the unit normal vector 
to the plane of the element. For 
a point on the exterior side of 
the element, we have the follow- 
ing signs on ZNP or Z for the 
primary and reflected elements: 

Sign of 

12 ZNP or Z 

1.3 + 

2.4 
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Square of the distance of the 
particle position from the ele- 
ment centroid. 

Square of the maximum diagonal of 
the element. 

Time step. 

ELEMENT NUMBER This is the sequence number of the 

element in the primary (i.e., unre- 
flected) surface description data 
set as input via the no. 12 cards. 

Calculation of the velocity continues, but return to the calling program 
is done with parameter INBODY set to unity instead of zero. 

If the calling program is SR TRAJEC, SR TRAJEC prints (in addition 
to and following the above): the particle coordinates inside the body at 

the end of the current time step, the initial y and z coordinates of the 
trajectory and the count, JT , of sequential trajectories that have re- 
sulted in impact. SR TRAJEC then calls SR IMPACT which may, at the 
discretion of the user, adjust initial trajectory coordinates,* and, if 
JT < JLIM (JLIM is set in a program above TRAJEC in the calling hierarchy), 
another trajectory is begun with the adjusted initial coordinates. 

Special Calculations 

Special calculation procedures must be used for points that are very 
close to element surfaces, to element edges, and for lifting elements, 
very close to extensions of element edges. Details are given in Appendix 
C. For nonlifting elements and for calculation of source contributions 
from lifting elements the special procedures have been found to work satis- 
factorily, and we expect no problems to be encountered by users. On the 
other hand, completely satisfactory and trouble-free procedures do not 
seem to be available for some dipole contribution calculations. For 
example, as a lifting element edge is approached during a particle trajectory 


ROSQ 

TSQ 

H 


*When used with PGM TANTRA, SR IMPACT must always be a dummy subroutine 
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calculation, anomalously large flow velocities may be returned by SR FLOVEL, 
which may cause the integration time step to fall to such a low value as to 
indicate a stall or stagnation condition. This usually can be taken to be 
tantamount to body penetration, however, the user would be well advised to 
examine the situation closely. In any case, when the integration time step 
falls below a threshold value set by SR T RAJ EC, a comment is printed and the 
trajectory is terminated. (See the discussion of SR TRAJEC , p. 74.) 

Array Dimensions 

The dimensions of arrays in COMMON/ BAS DAT/ must be adequate to accomo- 
date the data stored on unit 18 previously by program DUGLFT, and must be 
identical to those in the SR SETFLO COMMON/ BASDAT/. Details are included 
in the discussion of SR SETFLO below. COMMON/BASDAT/ appears in FLOVEL and 
SETFLO but nowhere else; neither subsidiary subroutines required for the 
FLOVEL calculations (VFNLFT, VFLIFT, NEARF AND WNEARF) nor the trajectory 
codes contain this labeled COMMON. 

SUBROUTINE SETFLO 

Subroutine SETFLO is called by programs FLOPNT , ARYTRJ , CONFAC and 
TANTRA prior to any calculation to recover data previously stored by DUGLFT 
on unit 18 (see pp. 55-57) that are needed by SR FLOVEL to calculate flow 
velocities about a three-dimensional body. To ensure that the proper set 
of data is recovered, SETFLO reads a four-character Hollerith body identifier, 
and checks to see if this is identical to the identifier obtained from the 
DUGLFT output (parameter CASE on DUGLFT input card no. 2). If not, a 
comment is printed and the calculation is stopped. 

COMMON/BASDAT/ contains all of the arrays that are recovered from the 
DUGLFT unit 18 output file. These arrays have variable dimensions in DUGLFT, 
and if dimensions are changed there, care must be taken to ensure that 
COMMON/BASDAT/ dimensions are adequate to contain the data. 
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Array dimensions in terms of the symbolic parameters used in DUGLFT 
are given in COMMENT statements at the beginning of the SR SETFLO FORTRAN 
listing. Definitions of the symbolic parameters are given above on pp. 43-44, 
and the contents of the unit 18 output file are described on pp. 55-57.. 

COMMON/BASDAT/ appears only in subroutines SETFLO and FLOVEL, but it 
must be identical in these subroutines. Therefore, if dimensions are 
changed in the SETFLO COMMON/BASDAT/ they must be changed identically in 
FLOVEL. 

If flow around the body is calculated by means other than the Hess 
method, SETFLO must be replaced by a dummy subroutine. 

PROGRAM FLOPNT 

General Description 

This program computes flow velocities at an array of points in three- 
dimensional space. The array is oriented parallel with the three coordinate 
axes discussed above on pp. 19-20. Flow velocities are computed by SR 
FLOVEL, which uses data that, for example, are prepared by program DUGLFT 
for flow about an arbitrary three-dimensional lifting body. 

Initial coordinates, array increment values along the three coordinate 
directions and the number of increments desired along each direction (in- 
cluding the initial point) are input. Also input are integers M(3) which 
control the order of incrementing along three axes. For example, suppose 
M(1 ) = 3, M(2) = 1, M(3) = 2: 

1. The x and z coordinates are held fixed while y is incremented over 
its range. 

2. y is returned to its initial value, z is incremented once, and y 
is incremented over its range. 

3. This is repeated until z covers its complete range. 
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4. z is returned to its initial value, x is incremented once, and y 
is incremented over its complete range. 

5. etc. 

The printed output is self-explanatory and consists of point coordinates, 
velocity components and speed. 

If data prepared by DUGLFT are used, SR SETFLO reads these data from unit 
18; units 5 and 6 are used for input and printing, respectively. 

Subroutines called are: SETFLO, FLOVEL. 
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FLOPNT Card Input 


Card 

No. Variables and Format Description 


1 

2 

3 

4 


5 

6 

3' 


KASE, (A4) Body identification. Read by SR SETFLO 

and must be identical to identifier on 
card 2 of DUGLFT. 


H0LL(18), (18A4) Run identification. 

M(3), (312) Coordinate incrementation sequence control. 

(See discussion above.) 


X( I ) , D( I ) , N( I ) : I = 1 

x(D 

(2E10.0, 14) 

(Col. 1-10) 


D(l) 


(Col. 11-20) 


NO) 


(Col. 21-24) 

X( I ) 5 D( I ) , N( I) : I = 2 

Same as card 

X( I ) , D( I ) , N( I ) : I = 3 

Same as card 


initial x coordinate \ 

Mdimensionless) 

x coordinate increment' 

number of increments desired in 
the x direction (including 
initial value). 

4 but for the y axis. 

4 but for fhe z axis. 


Cards 3-6 are repeated 
for another array. 


3 


Blank card 


A blank card 3 terminates the run. 
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TRAJECTORY CODE DESCRIPTIONS 


GENERAL UTILITY CODES 
Subroutine PARTCL 

Subroutine PARTCL is called by all three of the executive trajectory 
codes (Table IB) to input particle specification data and compute still-air, 
terminal particle settling speed and other data that depend on particle type. 
This is a particle type - specific code, the version used here being for 
water drops. It calls SR FALWAT. 


Subroutine TRAJEC 

Trajectories are calculated by sr TRAJEC with the assistance of: 

SR DVDQ, the numerical integrator code, SR FLOVEL and the functions PR FUN 
and IMPACT. It also stores trajectory point coordinates at user-specified 
(normalized) time intervals (TPRINT) in arrays XPL0T(60), YPL0T(60), 
ZPL0T(60), providing logical parameter IPLOT is specified as true. 

SR DVDQ uses a variable time step in integrating the particle equations 
of motion. An initial minimum time step, HMINI , is input by the user or set 
to 0.005 on default of input. (See program ARYTRJ input card 4.) During 
the calculations, the current value of the minimum time step is called 
HMIN and the time step is H. When DVDQ finds that H < HMIN, HMIN is set 
to H, a comment to this effect is printed, and the calculation is continued 
unless HMIN is found to have fallen to or below a threshold, HEPS, where 
HEPS = 2 x 10 5 HMINI. In this case, a comment, 

HMIN .LT. X.XXXXX E-XX THIS INDICATES STAGNATION 

is printed, the trajectory is terminated and control is returned to the 
calling program. 
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Function PRFUN 


Function PRFUN is a particle type - specific code which is called by 
T RAJ EC to provide the N D - N R relation used in calculating the particle 

equations of motion (eq. (1)). Actually, through use of the pre-calculated 

quantity COF (= ^Q jS v s ^p/^p s )> PRFUN returns the factor on the first term 

on the right side of eq. (1) which when multiplied by v - v yields the 

3 P 

particle equation of motion. The version of PRFUN used here is for water 
drops, and it calls functions CDRR and WCDRR. 


Subroutine IMPACT 

Subroutine IMPACT is called by TRAJEC following penetration of a 
particle into the body. When used with CONFAC, IMPACT is programmed by 
the user to adjust trajectory initial y and z coordinates such as to avoid 
impaction by the next trajectory, accordingly IMPACT is a problem-specific 
code. No such adjustment is required for cases run under control of 
ARYTRJ and TANTRA*, so that a dummy version of IMPACT is used. Examples 
of IMPACT use are as follows, in which we assume the coordinate system 
shown in Figure 2 and described on pp. 19-20. 

Suppose we wish to compute a trajectory to a target point on the 
starboard side of an airplane fuselage. Since all initial points of our 
trajectories will lie in the constant plane x = x^ (which is far ahead 
of the airplane in the unperturbed free stream), we need be concerned only 
with the initial y and z coordinates, y^ and z^ . If impaction on the body 
occurs, the most straightforward way to avoid impaction during the next 
trajectory is to adjust its y^ value to be further outboard in the star- 
board direction (see Fig. 2). Thus, to initiate the next trajectory we 
use y. . = y. + e , where e is a (normalized) positive increment of our 
choice. For example if e = 0.01, SR IMPACT would be: 


*Be very sure that IMPACT does not adjust initial trajectory coordinates 
during tangent trajectory determination under control of TANTRA. 
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SUBROUTINE IMPACT(YI ,ZI) 

YI + YI + .01 

RETURN 

END 

If our target point were to be near the underside of a wing, adjust- 
ment of initial coordinates would most likely be z- = z i - e, and if we 
choose e = 0.05, then SR IMPACT would be: 

SUBROUTINE IMPACT(YI.ZI) 

ZI = ZI - .05 

RETURN 

END 

Subroutine DVDQ 

This is the variable order, ordinary differential equation integrator 
of Krogh (ref. 7). Operating instructions, which have proven to be quite 
adequate, are found in the glossary of the DVDQ FORTRAN listing. This 
version automatically adapts to the word size of the computer used. 

A trajectory termination flag is set in SR DVDQ when the particle 
reaches a user-set x coordinate limit. This condition is found by interpola- 
tion (NGE = 0) or extrapolation (NGE =1). Interpolation is specified in 
the code as supplied, and is the more efficient and more accurate procedure. 
However, when termination is desired immediately upstream of the body sur- 
face, for example at a "leaking" aperture orifice (which would be paneled), 
it is necessary to use extrapolation to avoid overshoot and penetration of 
the body. Parameter NGE is specified in a DATA statement in SR TRAJEC. 

PROGRAM ARYTRJ 
General Description 

SR ARYTRJ is called to compute particle trajectories initiated at an 
array of points in three-dimensional space. Particle properties are com- 
puted by SR PARTCL and SR PRFUN. Flow velocities are computed by SR FL0VEL, 
which uses data that, for example, are prepared by program DUGLFT for flow 
around an arbitrary three-dimensional body. SR DVDQ integrates the particle 
equations of motion. 
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Initial coordinates, in the coordinate system defined above on pp. 19-20, 
of the initial point array, array increment values for the three coordinate 
directions and the number of increments desired along each direction (in- 
cluding the initial point) are input. Also input are integers M(3) which 
control the order of incrementing along the three axes and a skip parame- 
ter NSKIP. For example, suppose M(l) = 3, M(2) = 1, M(3) = 2: 

1. The x and z coordinates are held fixed while y is incremented over 
its range. 

2. y is returned to its initial value, z is incremented once, and y 
is incremented over its range. 

3. This is repeated until z covers its complete range. 

4. z is returned to its initial value, x is incremented once, and y 
is incremented over its complete range. 

5. etc. 

Trajectories are computed to the limiting x coordinate value XLIMIT or 
until penetration of the body is sensed. 

If not every trajectory is desired, the parameter NSKIP is set greater 
than zero. Then, after the first trajectory, only every NSKIP + 1 th tra- 
jectory is computed. 

Subroutines Required 

FLOVEL, SETFLO, PARTCL, FALWAT, T RAJ EC, IMPACT (dummy), PRFUN, 

DVDQ, WCDRR, CDRR 


External Storage Units 

Units 5 and 6 are the system input and print units, respectively. 

Unit 9 is used for temporary storage. 

Unit 10 is used to store trajectory data for plotting by PGM STEREO. 
Unit 18 is used by SR SETFLO for input of data prepared by PGM DUGLFT. 
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Printed Output 


The printed output is largely self-explanatory. For each trajectory 
are printed at time interval TPRINT : time, point coordinates (X, Y, Z), 

particle velocity components (VPX, VPY, VPZ), flow velocity components 
(VX, VY, VZ), time step (H), Reynolds number (R), acceleration modulus (AC) 
and cumulative number of flow velocity computations (NEVAL). (All dimen- 
sionless) 

Other quantities are: angle between the projection of the initial 

flow velocity vector on the z = 0 plane and the x axis (ALPHAO), angle 
between the initial flow velocity vector and its projection on the z = 0 
plane (BETAO) , angle between the projection of the final particle velocity 
vector on the z = 0 plane and the x axis (ALPHAR), angle between the final 
particle velocity vector and its projection on the z = 0 plane (BETAR), 
direction cosines of the drag vector at the final point, and the angle 
between the projection of the drag vector on the z = 0 plane and the x 
axis (A), and the angle between the drag vector and the z axis (GAMMA). 

(All angles are in degrees.) 
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ARYTRJ Card Input 


Card 

No. Variables and Format 

1 KASE, (A4) 


2 HOLL (18), I PLOT, 
( 1 8A4 , 7X, LI) 


3 V, ELL, RHO, TEMP, 
XFINAL, (8F10.5) 


4 TPRINT , HI, HMINI, 
EPS 1(3) , (8F10.5) 


5 DIAM, ( FI O.O) 

6 M( 3) , NSKIP, (414) 


7 X( I ) , D(I), N( I ) ; 

I = 1 (2F10.0, 14) 


8 X( I ) , D(I), N(I) ; 
I = 2 

9 X( I ) , D( I ) , N(I) ; 
I = 3. 


Description 

Body identification. Read by SR SETFLO. Must be identical 
to parameter CASE on card 2 of the DUGLFT input. 


HOLL 72 columns of Hollerith run identification. 

IPLOT Logical variable: if true, trajectory data 

(col. 80) are written on unit 10 for plotting by PGM 
STEREO. 


V 

(col. 1-10) 
ELL 

(col. 11-20) 
RHO 

(col. 21-30) 
TEMP 

(col. 31-40) 

XFINAL 
(col. 41-50) 


Free stream airspeed (m s” 1 ) 

Characteristic dimension of the body (m). 
Corresponds to L as defined for eq. (1). 

Ambient air density (kg m” 3 ) 

Ambient temperature (°K) 

x coordinate for trajectory cut off (i.e., 
maximum x coordinate) (normal ized , dimensionless) 


TPRINT Time interval for trajectory point print. 

(col. 1-10) Default value = 0.1. 

HI Initial numerical integration time step. (See 

(col. 11-20) SR DVDQ.) Default value = 0.1 

HMINI Initial numerical integration minimum time 

(col. 21-30) step. (See SR DVDQ.) Default value = .005. 

EPSI (3) Parameters used to control numerical integra- 
ted. 31-60) tion local error. (See SR DVDQ.) Default 

values = 1 .OE-5. 

All normalized dimensionless. 


Water drop diameter (ym). This card is read by SR PARTCL. 

M ( 3 ) Array incrementation control. 

(col. 1-12) 


NSKIP 

(col. 13-16) 

X(l) 

(col. 1-10) 

D(l) 

(col. 11-20) 

NO) 


Array skip parameter. (See discussion above.) 


Initial x coordinate 


x coordinate increment 


(dimensionless) 


Number of increments desired in the x direction 


(col. 21-24) (including the initial value). 
Same as card 7, but for the y direction. 


Same as card 7, but for the z direction 


5' Cards 5 - 9 are re- 
peated for another 
particle and another 
array 


5 Blank card A blank card 5 terminates the run. 
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PROGRAM CON FAC 


General Discussion 

Program CONFAC computes trajectories to user-specified target points. 
It operates in two modes: 


1. Single trajectories are calculated to each target point (NW = 0). 

2. A central trajectory is computed to the target point, and NW 
trajectories, evenly spaced about a circle in the target plane 
of radius RW about the central trajectory, are calculated such 
as to define a particle flux tube. 


Mode 2 is used to calculate concentration factor, Cp, which is the ratio of 
particle flux at the target point to the free-stream particle flux. It is 
easy to show that 


C 


F 


area of flux tube cross section in the free stream 
area of flux tube cross section at the target point 


The areas are those of plane polygons of NW vertices as calculated by 
SR POLYGO. Concentration ratio, C M , the ratio of particle concentra- 
tion at the target point to free stream concentration, is obtained via 
the relation 


C 


M 



The desired trajectories are calculated by an iterative method which 
finds a trajectory that passes within a user-specified distance-tolerance 
(RW*T0L) of the desired target point. To initialize, the user may input 
four sets of coordinate guesses: two sets of y and z coordinates for the 

initial and target planes. These two sets of data are requried by the 
iterative trajectory method since it always must have at least two sets 
of initial and final y and z coordinates in storage. No special care need 
be taken in making these guesses since convergence should be rapid as long 
as the coordinates are in the correct general neighborhood. On default of 
input, these coordinate guesses are supplied by the code, and these default 
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values should be adequate. However, should the user wish to restart a 
trajectory iteration, he can initialize with data sets from the output of 
the previous calculation or any values of his choice, such as to avoid an 
initiation that is identical to the previous one. 

The trajectory iteration procedure is described in detail in refer- 
ence 16. (See pp. 13-16 and Appendix A of ref. 16.) SR MAP controls 
the iteration and calls SR TRAJEC to calculate trajectories. If convergence 
is not achieved after calculating twenty-five trajectories, the calcula- 
tion proceeds to the next particle or stops. The limiting number of tra- 
jectories can be changed by changing the value of HIM in a DATA statement 
in SR MAP. 

SR IMPACT is a problem-specific code whose purpose is to adjust tra- 
jectory initial y and z coordinates when penetration of the body occurs such 
that penetration will be avoided on the next attempt. After twenty-five 
penetrations, the calculation proceeds to the next particle or stops. 

The limiting number of penetrations can be changed by changing the value 
of JLIM in PGM CON FAC. 


Subroutines Required 

FLOVEL, MAP, PARTCL, POLYGO, DVDQ, SETFLO, FALWAT , PRFUN, IMPACT, 
TRAJEC, TRANSF, MATINV , WCDRR, CDRR. 


External Storage Units 
Same as for ARYTRJ. 


Printed Output 

The printed output is largely self-explanatory, and contains all of the 
data described for PGM ARYTRJ. 
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Detailed trajectory data are printed only for final trajectories which 
result from successful convergence to desired target points. All coor- 
dinates and velocities in this printout are defined in terms of the "flow 
coordinate system", which is as described above on pp. 19-20. 

For trial trajectories, only initial and final y and z coordinates 
are printed. Except for the initial coordinate guesses, these coordinates 
are given in the "flux tube coordinate system", and are so labeled in the 
printout. The flux tube coordinate system is different from the flow sys- 
tem only for flux tube peripheral trajectories. Then at every point along 
the flux tube, the y-z plane of the flux system is normal to the central 
trajectory, and its origin in the y-z plane is at the central trajectory. 

In cases where a flux tube central trajectory is not defined (i.e., for 
single trajectories (NW = 0) and during calculation of a central trajectory) 
the flux tube and flow coordinate systems are identical. 

When flux tubes are calculated, a summary of results is printed. This 
includes : 

1. Initial and final coordinates of the central trajectory, de- 
fined in the flow system. 

2. Coordinates of the polygon vertices (i.e., peripheral trajectory 
intersections) in the initial and target planes, defined in the 
flux tube system. 

3. Polygon areas in the initial and target planes. 

4. Concentration factor and concentration ratio. 

Also printed are angles a.., 3 . (ALPHAO, BETAO) at the initial plane 
and a t , 3 t (ALPHAR, BETAR) at the target plane, which can be used to trans- 
form coordinates between the flow and flux tube systems. Given coordinates 
x, y, z in the flow system along with central trajectory initial (or final) 
coordinates x^ , y i , z^ (or x t , y t> z^) also in the flow system, we define 
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Ax = x - x i (or Ax = x - x t ). Ay = y - y. (or Ay = y - y t ) and Az = z - z . 

(or ha. - z - z^) , and designate x^, y^, z^ to be the corresponding coordinates 
in the flux tube system. Then 

x f = Ax cosa cose + Ay sina cose + Az sine 

y^. = -Ax sina + Ay cosa 

z f = -Ax cosa sine - Ay sina sine + Az cose 


and 


ax = x^ cosa cose - y^ sina - z^ cosa sine 

Ay = x^r sina cose + y^ cosa - z^ sina sine 

Az = x f sine + z f cose 

Here a is the angle between the projection of the velocity vector in the 
x-y plane and the x axis, and e is the angle between the velocity vector and 

its projection in the x-y plane, where the x-y plane is defined in the 

flow coordinate system. 
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CON FAC Card Input 



Card 

No. 

Variables and Format 


Description 

1 

KASE, (A4) 

Body identification. Read by SR SETFLO. Must be identical 
to parameter CASE on card 2 of the DUGLFT input. 

2 

HOLL ( 18) , I PLOT, 

HOLL 

72 columns of Hollerith run identification. 


( 1 8A4 , 7X, LI) 

IPLOT 
(col. 80) 

Logical variable: if true, trajectory data 
are written on unit 10 for plotting by PGM 
STEREO. 

3 

V, ELL, RHO, TEMP, 
XSTART, (8F10.5) 

V 

(Col. 1-10) 

Free stream airspeed (ms 1 ) 



ELL 

(Col. 11-20) 

Characteristic dimension of the body (m). 
Corresponds to L as defined in eq. (1). 

It is used to normalize all lengths. 



RHO 

(Col. 21-30) 

Ambient air density (kg m” ). 



TEMP 

(Col. 31-40) 

Ambient temperature (°K). 



XSTART 
(Col. 41-50) 

Initial x coordinate of trajectory, (normalized 
dimensionless) 

4 

TPRINT, HI, HMINI , 
EPSI (3) , (8F10.5) 

Same as for 

ARYTRJ . 

5 

NW, RW, TOL, (110, 
7FI0.5) 

NW 

(Col. 1-10) 

Number of trajectories used to define flux 
tube peripheries for concentration factor 
calculation. If NW = 0, single trajectories 
are calculated to target points defined by 
cards 9.* Right justify in the field. 



RW 

(Col. 11-20) 

Flux tube radius (normalized). 



TOL 

(Col. 21-30) 

Tolerance factor for trajectory iteration 
convergence. Convergence occurs when a 
trajectory passes within a distance of 
RW*T0L of the target point. 

6 

YE, ZE, YI, ZI 
(8F10.5) 

YE, ZE 

Initial guesses of trajectory y and z 
coordinates in the target plane 

7 

YE, ZE, YI, ZI, 
(8F10.5) 

YI, ZI 

Initial guesses of trajectory coordinates 
in the initial plane 



These coordinates are in the coordinate system used to 
define the body and the flow field (normalized, dimension- 
less). The data in the two cards can be very approximate, 
but if not blank, the two cards should not be identical. 

On input of two blank cards, the code supplies default 
estimates based on the first target coordinates. 

8 

DIAM, (7F10.0) 

Water drop diameter (pm). This card is input by SR PARTCL. 

^Parameters RW and TOL must be 

given non-zero values even when NW = 0. 
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CONFAC Card Input (cont.) 


Card 

No. Variables and Format __ Description 


9 


8 ' 


9' 


XW, YW, ZW, (8F10.5) 

Cards 8 and 9 are re- 
peated for as many 
particles as desired.* 


X I y I Z . C !°[ d i nateS ° f the target P° int (normalized, dimensionless). 


8 Blank card 


A blank card 8 terminates the run. 


*Previous trajectory y and z coordinates are used as trajectory iteration initialization estimates for 
each new target point. Thus, if target points are widely spaced, separate runs should be made for each. 
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PROGRAM TANTRA 
General Discussion 

The purpose of this code is to compute tangent particle trajectories to 
a three-dimensional body. The code is designed to be as general and as auto- 
matic as practical, but owing to the unlimited number of geometrical possi- 
bilities in three dimensions, some compromise is necessary. Since we cannot 
know a priori what parts of a body the tangents will touch, we do not, in 
general, have the option of specifying target points on the body or even 
target planes through the body. Therefore, we specify curves in the free 
stream well ahead of the body on which all trajectories are initiated for a 
particular tangent determination. 

Given the equation of the starting-point curve and an initial point on 
the curve, the code computes the trajectory from this point toward the body 
until penetration of the body occurs or until a specified x-coordinate stop 
point is reached. If penetration occurs, a specified coarse step is taken 
along the starting-point curve such that the resulting trajectory will pass 
further from the body, and another trajectory is computed. If penetration 
does not occur, the coarse step is taken along the starting-point curve in 
direction such that the resulting trajectory will pass closer to the body, 
and another trajectory is calculated. Once penetration occurs for trajec- 
tories that initially miss the body, or the reverse for trajectories that 
initially impact, the initial point is backed up one step along the start- 
ing-point curve, and the process of stepping toward or away from the 
body is resumed with a fine step size until the tangent trajectory is 
found. Thus, the tangent trajectory misses the body by no greater than the 
tolerance implied by the fine step size. Note that this does not imply 
that the tolerance is the fine step size. Separation of trajectories in 
the free stream will not be the same as separation of the same trajectories 
near the body, nor even approximately the same except for very large, heavy 
particles which have sufficient inertia to essentially ignore the flow 
around the body. In general, trajectory separations near the body will be 
less than in the free stream. 
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Specification of the starting-point curve is done via SR STRPNT, which 
in the version supplied uses straight line curves. The user provides the 
coordinates of two points on the line: 

Point 1. Initial coordinates of the initial trajectory. 

Point 2. Coordinates of any other point on the line chosen such 
that a small step taken along the line from point 1 
toward point 2 will initiate a trajectory that will 
pass closer to the body surface, or intersect it at 
a point that will project closer to the body center, than 
the trajectory initiated at point 1. 

This definition of Point 2 relative to Point 1 is such as to ensure that 
stepping along the starting-point curve will proceed in the proper direc- 
tion. Point 2 need not be, and in general will not be, the initial 
point of an actual trajectory. Note that both of these points must be 
sufficiently far upstream to be essentially in the free stream. Also 
specified are the coarse and fine stepping distances. All coordinates 
and distances are normalized. (See eq. (1).) 

If so specified (IPLOT = true), tangent trajectory data are stored 
on unit 10 for plotting later by PGM STEREO. 


Subroutines Required 

FLOVEL, SETFLO, PARTCL, FALWAT, STRPNT, T RAJ EC, IMPACT*, PRFUN, 
DVDQ, WCDRR, CDRR. 

External Storage Units 
Same as for ARYTRJ. 


*Be sure that IMPACT is a dummy subroutine. Resetting of initial trajectory 
coordinates by SR IMPACT will ruin a tangent trajectory determination. 
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Printed Output 


Trajectory data are as described for PGM ARYTRJ and are printed for 
every trajectory computed regardless of whether or not the trajectory is 
accepted as the tangent trajectory. Beyond that, the output is fully 
labeled and self-explanatory. All input data are printed: including the 

points used to define the starting-point line, coarse and fine increments, 
and starting point coordinates. The switching from coarse to fine step 
size is clearly identified, as are the tangent trajectory data. 
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TANTRA Card Input 


Card 

No. Variables and Format Data Description 


1-4 


Cards 1 through 4 are the same as for 
ARYTRJ. 


5 

6 


7 


DIAM, (7F10.0) 

DCORS, DFINE , (8F10.0) 


X,Y,Z,X1,Y1,Z1 , 
(8F10.0) 


Water drop diameter (pm). This card 
is read by SR PARTCL. 

Respectively, the coarse and fine step 
sizes to be used in stepping along the 
starting-point line (normalized, dimen- 
sionless). Card 6 is read by SR STRPNT. 

X,Y,Z Coordinates of Point 1, which 
specifies the initial trajec- 
tory coordinates on the start- 
ing-point line (normalized, 
dimensionless). 

XI ,Y1 ,Z1 Coordinates of Point 2, which 
is any point on the line that 
will initiate a trajectory 
that will pass closer to the 
body than the trajectory initia- 
ted by Point 1. (normalized, 
dimensionless) . 

Note that the starting-point line should 
be far enough upstream of the body to be 
essentially in the free stream. Card 7 
is read by SR STRPNT. 


6' Cards 6 and 7 are repeated for as many 

7' trajectories as desired. 


6 Blank card A blank card 6 signals end of calculation 

for this water drop, and another card 5 
is read. 


5' 

6 " 
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Tantra Card Input (cont.) 


Card 

No. Variables and Format Data Description 


7" 


6 Blank card 


5 Blank card A blank card 5 terminates the run. 
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PROGRAM STEREO 


General Discussion 

Program STEREO is used to plot results of the trajectory calculations. 
Both body and trajectories are plotted. The body data are obtained from 
unit 9, on which the data were stored by SR PINPUT under control of PGM 
PBOXC, and the trajectory data are obtained from unit 10, on which the data 
were stored under control of either ARYTRO, CONFAC or TANTRA. 

Plots are prepared in pairs, members of a pair being separated by a 
specified angle on each side of a specified viewing direction. Proper 
specification of the angles, which usually requires some trial-and-error- 
experimentation, may provide plots which can be used for stereo viewing as 
illustrated by Figure 9. 

The viewing direction is defined by specifying two angles, THETA and 
PSI. The operation of these angles is as follows: We assume a right-handed 

coordinate system with its positive z axis directed upward and the free- 
stream flow in the direction of the positive x axis. First rotate the 
coordinate system about the y axis by angle THETA such that positive THETA 
tilts the positive x axis upward. Then rotate about the new z axis by angle 
PSI such that for positive PSI the rotation is clockwise when viewed from 
above. The view direction separation angle, DELTA is applied to angle PSI 
such that the members of a stereo pair are actually viewed from angles 
THETA, PSI-DELTA and THETA, PSI + DELTA, and are plotted in that order. 

For a particular case (i.e., body and set of trajectories), the user 
must specify the number of trajectories and the (upstream) x coordinate at 
which plotting of the trajectory data is to be begun. This need not have the 
same value as the initial x coordinates of the data stored on unit 10. 

Translating and scaling of the data such that it will properly fit into 
the plot area is handled automatically by the program. 

Only system and plot subroutines are required. 
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Figure 9. Stereographic plot of an eight-trajectory, 20 ym-diameter water drop flux tube 
to a particle replicator mounted on the forward fuselage of a Lockheed C130A 
airplane. The central trajectory is also shown. C = 1.15. (Ref. 16) Three 

dimensional perspective can be attained by staring at the center of the figure 
and then crossing the eyes such that the two images merge. 


External Storage Units 


Units 5 and 6 are the system input and print units. 

Unit 9 contains the three-dimensional body surface data, plus some 
scaling information, as stored by SR PINPUT. 

Unit 10 contains the trajectory data as stored under control of PGMS 
ARYTRJ , CONFAC or TANTRA. 


Printed Output 

The printed output is simple. It consists of a run identification, the 
input data and some scaling information. For each trajectory is printed: 

1. the coordinates (XTRAJ, YTRAJ, ZTRAJ) of each point before translation, 
scaling and projection onto the plot plane, and 2. the translated, scaled 
and projected coordinates (XPLOT, YPLOT) of each point plotted. 



STEREO Card Input 


Card 

No. Variables and Format Description 

1 H0LL(18), (18A4) 72 columns of Hollerith run identification. 


2 


ICRT, NTRJS, X START, 
(LI, 19, F10.0) 

ICRT 
(Col .1 ) 

A logical variable which when 
true causes plotting to be 
via CRT. Otherwise, plotting 
is via pen and ink. 


NTRJS 
(Col. 2-10) 

Number of trajectories to be 
plotted. 


XSTART 
(Col. 11-20) 

x coordinate at which trajec- 
tory plotting is to begin. 

Need not correspond to the 
initial x coordinates of 
trajectories stored on unit 10. 


3,4 LINE1, LINE2, 
(7A6/7A6) 


Cards 3 and 4 are read only if ICRT is 
true (card 2). Two lines of 42 columns 
each of Hollerith labeling for a micro- 
fiche film. 


5,6 THETA, PSI, DELTA, 

HLABEL(18), (3F10.2/ 
18A4) 


THETA Viewing angles and viewing angle 

PSI separation (degrees). (See 

DELTA definitions above. ) 

HLABEL 72 columns of Hollerith labeling 

for the plots. 


5', 6' Cards 5 and 6 are repeated for as many additional plot 
pairs as desired. 


5,6 Blank cards 


Blank cards 5 and 6 terminate the run. 
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VALIDATION 


PRIOR WORK 

Hess and Smith report a large number of validation studies for non- 
lifting flow (refs. 4 and 5). These studies include flow about bodies 
of simple geometries for which analytical solutions are known, and 
comparison of calculated results with wind tunnel data for complex air- 
craft bodies. Excellent agreement is obtained in almost all cases. 

Hess (ref. 3) presents the results of a variety of studies on 
airfoils and combinations of airfoils with nonlifting bodies using his 
lifting code. Pressure distributions, both chordwise and spanwise, are 
generally in good agreement with measured data. Local lift coefficient, 
however, is always overpredicted, though there is usually good agreement 
in shape of the calculated and observed distribution curves. Over- 
prediction of lift coefficient is substantial, being typically in the range 
of 10-20%. Results of one study in which boundary layer effects were 
approximately accounted for indicate that most, if not all, of this cal- 
culation error is caused by neglect of viscosity (ref. 3). In some cases, 
neglect of local compressibility also substantially affects the lift cal- 
culations, but Hess concludes that a simple correction suitable for lifting 
flow is not available (ref. 3). 

Studies of trajectory calculations using the nonlifting code are 
discussed in references 2 and 16, and those results are merely summarized 
here: (1) Very small particles were calculated to follow flow streamlines 

around an ellipsoid with negligible deviation, and from this we conclude 
that our numerical integration procedure is highly accurate. (2) Trajectory 
results for water drops over a range of sizes in flow about an analytical 
ellipsoid were compared with results from an ellipsoid defined by quadri- 
lateral elements. Small differences were found except for a drop on the 
edge of a shaded zone, where numerical difficulties are inevitable in any 
case. (3) Tangent trajectory results calculated for flow about ellipsoids 
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were compared with prior calculated results, and acceptable agreement was 
found. (4) Tangent trajectory results were compared with observed data 
and prior calculations for four wind tunnel tests of flow about an ellip- 
soidal forebody. Maximum impact distance along the surface from the nose 
point was consistently overpredicted by about 60%, but this is an improve- 
ment over the prior calculations which overpredicted in the range from 
60% to 160%. Moreover, overprediction by the prior study consistently 
increases as droplet size decreases, whereas this trend is not particularly 
evident in our calculations. Since the trajectory calculations become 
more demanding as particle size decreases, this improvement is suggestive 
of improved capability. 


ADDITIONAL VERIFICATION 

NACA TN 3839 (ref. 17) reports wind tunnel data for impact of water 
drops on several airfoils. Calculated results also are given for the 
NACA 65-j -212 airfoil at 4° angle of attack, so that we have selected this 
case for verification study, plus the same airfoil at 0° angle of attack. 
Airfoil coordinate points (Table 4a) were taken from Table 1 of NACA TN 
3047 (ref. 18), and plots of the digital representation of the airfoil are 
given in Figure 10. Also in TN 3839 are wind tunnel data for the NACA 
65 i- 212 airfoil that was yawed at an angle of 35° to the onset flow such 
as to represent an airfoil with a 35° sweep angle. The airfoil crossection 
in the plane normal to its leading edge is the same as for the unswept 
airfoil. To determine the crossection coordinates in the plane parallel 
with the onset flow, simply multiply the z coordinates in Table 3a by 
cos(35°), and these coordinates are listed in Table 4b. Figure 11 shows 
a plot of the swept wing digital description. 
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TABLE 4 


a. NACA 65i-212 AIRFOIL CROSSECTION COORDINATE POINTS (Ref. 18) 



Lower Surface 



Upper Surface 


X 

z 

s* 

X 

z 

s* 

0 

0 

0 

0 

0 

0 

.00577 

-.00870 

.0100 

.00423 

.00970 

-.0100 

.02609 

-.01686 

.0325 

.02391 

.02058 

-.0325 

.07627 

-.02745 

.0830 

.07373 

.03593 

-.0840 

.15121 

-.03727 

.1600 

.14879 

.05073 

-.1610 

.25094 

-.04510 

.2590 

.24906 

.06300 

-.2620 

.35058 

-.04882 

.3590 

.34942 

.06942 

-.3630 

.45019 

-.04854 

.4605 

.44981 

. 07044 

-.4640 

.54983 

-.04317 

.5610 

.55017 

.06507 

-.5705 

.64957 

-.03351 

.6625 

.65043 

.05411 

-.6720 

. 74947 

-.02164 

.7640 

.75053 

.03954 

-.7730 

.84955 

-.00956 

.8640 

.85045 

.02302 

-.8750 

.94983 

-.00039 

.9650 

.95017 

.00671 

-.9775 

1 . 00000 

0 

1.0150 

1.00000 

0 

-1.0275 


*s is the distance from the leading edge along the airfoil surface. 


b. CROSSECTION COORDINATE POINTS FOR NACA 65 j -212 
AIRFOIL WITH 35° SWEEP 


Lower Surface Upper Surface 


X 

z 

X 

z 

0 

0 

0 

0 

.00577 

-.00713 

.00423 

.00795 

.02609 

-.01381 

.02391 

. 01 686 

.07627 

-.02249 

.07373 

.02943 

.15121 

-.03053 

.14879 

.04156 

.25094 

-.03694 

.24906 

.05161 

.35058 

-.03999 

.34942 

.05687 

.45019 

-.03976 

.44981 

.05770 

.54983 

-.03536 

.55017 

.05330 

.64957 

-.02745 

.65043 

.04432 

.74947 

-.01773 

.75053 

.03239 

.94955 

-.00783 

.85045 

.01886 

.94983 

-.00032 

.95017 

.00550 

1.00000 

0 

1.00000 

0 



ing edge 




Upper surface 


m 




^Trailing edge 


Lower surface 


Figure 11. Computer-drawn plots of the digital 
description of the NACA 65 x -212 air- 
foil with 35° sweep angle. 


’VTrailing edge 



Data were collected in the Lewis Laboratory Icing Research Tunnel for 
the following atmospheric conditions (ref. 17): 

Airspeed 152 kts (78.195 ms 1 ) 

Temperature 50°F (283.16°K) 

Pressure 28.1"Hg (94985 Pa) 

Air density and viscosity are calculated to be 1.1686 kg m and 1.76520X10 
kg m" 1 s" 1 . 

As shown in Figure 10, the unswept airfoil is represented by twenty 
lifting strips, each with twenty-six on-body elements and three wake 
elements, including the semi-infinite final wake element. A (plus) y=0 
symmetry plane was defined which has the effect of doubling the airfoil 
to forty lifting strips. The swept airfoil (Figure 11) is represented 
by nineteen lifting strips, each with twenty-six on-body elements and 
three wake elements, including the semi-infinite final wake element. No 
symmetry is used. The pressure equality Kutta condition and the step 
function spanwise vorticity variation options were exercised for 
both airfoils. Chord lengths were taken to be 13 and 87.5 inches re- 
spectively for the unswept and swept airfoils. 

Data were collected in the wind tunnel using water sprays with three 
different distributions of water droplet sizes (ref. 17). To compare with 
tangent trajectory calculation results, we are interested in the maximum 
distance aftward from the leading edge, on both the upper and lower surfaces, 
at which droplet impact was observed. Such distances measured along the 
surface of the airfoil, s and s , correspond approximately to impact points 

U X- 

of barely subtangent trajectories of droplets of maximum diameter, S max , in 
the sprays. For the unswept airfoil, we calculated s u and for each of 
the three sprays for 0° and 4° angles of attack. For the swept airfoil, we 
calculated s^ and s^ for the spray of intermediate droplet size for angles 
of attack of 0° and 4.3°. 

Results are listed in Table 5 and displayed graphically in Figure 12 
for the unswept airfoil. Table 6 lists results for the swept airfoil. 
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TABLE 5 

TANGENT TRAJECTORY RESULTS FOR THE UNSWEPT AIRFOIL 
a. 0° Angle of Attack 


6 max 

(ym) 

Observed 

— s /s ► 

u z New 

Calculation 

% Error 

59 

.250/. 202 

.230/, 260 

-8/29 

48 

.199/. 170 

.240/. 185 

21/9 

29 

.121/. 081 

.144/. 096 

19/19 




cr 

-C* 

o 

Angle of Attack 


r 


— V s * 

% 

Error 

max 


New 

NACA New 

NACA 

L*!lL 

Observed 

Calculation 

Calculation Calculation 

Calculation 

59 

.109/. 460 

.160/. 456 

.20/. 45 47/-0.9 

83/-2 

48 

.085/. 394 

.082/. 452 

.17/. 42 -4/15 

100/7 

29 

.041/. 248 

.036/. 359 

.10/. 33 -12/45 

144/33 


Note: All s values are normalized by dividing by the chord length, 13". 


TABLE 6 

TANGENT TRAJECTORY RESULTS FOR THE AIRFOIL WITH 35° SWEEP 


max 


= 48 pm 


Angle of 
Attack 

Observed 

— V s *— 

Calculated 

% Error 

0° 

.05/. 04 


.073/. 030 

46/ -25 

4.3° 

.01/. 14 


.009*/. 156 

-10/11 


Note: All b values are normalized by dividing by the chord length, 87.9 !! . 


*The particle stalled at a distance of about .012" from an element edge. 

See the discussion of special calculations (pp. 69-70), p. 74 and Appendix C. 
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.3 


.2 


.1 



max. 


Figure 12. Calculated and observed results for water drop impaction on 

an unswept NACA 65i-212 airfoil. s u and s^ are maximum impact 
distances from the leading edge, measured along the airfoil 
surface, on the upper and lower sides of the airfoil respec- 
tively. <5 is the diameter of the largest drop in the water 



As might be expected, agreement with the observed data is not outstanding. 
However, when experimental problems which must have influenced the observed 
results are considered (speculatively), and also we consider that the 
use of rather coarsely resolved flat elements to represent the slightly 
curved airfoil surface will at least occasionally cause substantial error 
in determining impact points of nearly tangent trajectories, we conclude 
that the comparisons are satisfactory. Moreover, for the unswept airfoil 
at 4° angle of attack, we can compare the new calculation results with 
the calculated results reported in reference 17. The average absolute 
percent error for the new calculation is 21% while it is 47% for the 
NACA calculations. Thus, the new code appears to provide substantial im- 
provement over the previous calculations. 

CALCULATION TIMES 
IBM 370/3033 Computer 

Element processing (program DUGLFT) and tangent trajectory (sub- 
routine TANTRA) calculations for the unswept NACA 65 i- 212 airfoil were 
performed on the IBM 370/3033 computer at the NASA Lewis Research Center. 
The DUGLFT run to process the 1200 elements at 0° angle of attack re- 
quired 197 seconds of CPU time. Again for the 0° angle of attack case, 
tangent trajectories were computed in two runs: one for all three par- 

ticles to the upper airfoil surface, and the other for the three particles 
to the lower surface. Results are as follows: 

a. Upper Surface 

5 max Number of Number of 

ImlL Trajectories Velocity Calculations 


59 

10 

3411 

48 

13 

4123 

29 

13 

A n no 
*+/ C.C 

Total CPU time = 

2892 seconds. 
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b. Lower Surface 


6 max 

Number of 

Number of 

(ym) 

Trajectories 

Velocity Calculations 

59 

9 

2422 

48 

10 

3249 

29 

13 

5039 


Total CPU time = 2492 seconds. 

Each trajectory calculation is initiated at x = -5 and stopped at x = 1 
unless impaction occurs. (Distances are normalized by dividing by the 
chord length.) 

Since most of the computing time is used for velocity calculation, 
we have included the number of velocity calculations in the tabulation 
above, and a rough value of the computing time per velocity calculated 
from these data are (for 1200 lifting elements): 

upper surface 0.236 second/velocity 

lower surface 0.233 second/velocity. 

CPC 6600 Computer 

Calculations similar to those described above were done for the 
NACA 65i -212 airfoil with 35° sweep on the CDC 6600 computer at the Air 
Force Geophysics Laboratory, Hanscom AFB, Massachusetts. 

PB0XC calculations required to prepare plots of the digital descrip- 
tion of the airfoil, which is comprised of 551 elements, required 22 
seconds of CPU time. (PB0XC processing of the 1200 elements of the un- 
swept wing required 21 seconds of CPU time on the CDC 6600 computer.) 

Trajectory calculations, for the water drop of diameter 48 ym, 
required the following: 
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Angle of 
Attack 
(deg. ) 

Upper or 

LU'VYCI 

Surface 

MnmKov' n* f 

MMIIIk/OI V 1 

Trajectories 

Number of 
Velocity 
Calculations 

CPU 
Time 
(sec. ) 

0 

4.3 

0 

4.3 

Upper 

Upper 

Lower 

Lower 

5 
7 
7 

6 

2824 

3058 

3514 

4245 

1039 

1102 

1273 

1588 


From these data we compute the following rough estimates of comput- 
ing time per velocity (for 551 elements): 

upper surface 0.364 second/velocity 

lower surface 0.369 second/velocity 

Discussion 

As discussed on p. 13 above, the complexity of, and hence the time re- 
quired for, calculation of a velocity contribution from a particular element 
increases as separation between the calculation point and the element 
decreases, according to which of three classifications, far field, inter- 
mediate field or near field, the separation distance falls into. Since 
in the tangent trajectory calculations described above the distance of 
each element from each calculation point changes with each time step during 
integration of the particle equations of motion, we have a changing mixture 
of far field, intermediate field and near field velocity contribution calcula- 
tions, so that it would be difficult to extract from the results generally 
meaningful information on calculation time required per element per velocity. 
The reader should be aware that the calculation time is sensitive to the 
mixture of far, intermediate and near field elements (as well as total number 
of elements), and this mixture and its change with time during trajectory 
calculation will be determined by the geometry of the particular problem. 

Calculation time is also sensitive to particle size, particularly in 
curved flow fields. For large (i.e. , massive) particles, which possess 
sufficient inertia to substantially ignore complexities in the flow field, 
trajectory calculations require few velocity calculations and hence little 
computer time. On the other hand, very small particles closely follow the 
flow streamlines, which requires small time steps and much more computer 
time. 
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EXAMPLE PROBLEM 


GENERAL DISCUSSION 

Example card input data and printouts are given below for programs 
DUGLFT and FLOPNT. Example problem input data, printouts and computation 
times for the other codes are given in reference 2 for the nonlifting 
code. 

THE TEST BODY 

The test body is a simple, hollow, box fuselage with airfoils attached 
(Fig. 13). The body is constructed of two nonlifting sections (fuselage), 
and one lifting section (airfoil), with one symmetry plane. The airfoil 
consists of four lifting strips, including an extra strip that extends 
through the fuselage to the symmetry plane, each with ten on-body elements 
plus three wake elements. This is the same test body described by Mack 
for the McDonnell Douglas lifting code (ref. 6). 

EXAMPLE PROBLEM CALCULATION 

The example calculations are for the test body with a ten-degree nose- 
up angle of attack. The step function and pressure equality options are 
selected, and semi-infinite final wake elements are used. Flow velocities 
at six off-body points are calculated by PGM DUGLFT, and calculations for 
these same points are repeated by PGM FLOPNT. 

The data set stored on unit 18 by DUGLFT may be saved (catalogued) 
such as to be available to FLOPNT during a separate run. 

CALCULATION TIMES 

CPU times for the test problem are (seconds): 

CYBER 750 


PBOXC 

1.27 

DUGLFT 

4.48 

FLOPNT 

0.23 
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Figure 13. Test problem body for the Hess lifting 
code (ref. 6). 
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STORAGE 


High-speed access storage required fo«* the test problem are (decimal) 



CYBER 750 


(kilo words) 

PROXC 

27 

DUGLFT 

29 

FLOPNT 

18 
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EXAMPLE PROBLEM CARD INPUT 


Id 


i a! 5 £S^ 5 £ 5 £S^!£ 3 Ei 

I aE]IBaBaElgBgKIEIEKI 


lEBCBKBKBRBILil 

ICBEaCBCBUQI 

IBBBBBBBBBBBI 

IL’BFBfnBLnBFiBL'lj 

HUBEiBEtiBrcBrcBI»i 

I^BEBBfllrlBKfllCil 


IBBBBBBBBBBBI 

ibbbbbbbbbbbbE 




1151 
IKSl 
■ehB 
zM&M 

jSbbibi 

SBI3BBI 


I^BBBBIBBB 
■t^BBBBBBBB 
?BF1BBB ■ 

SBBBBBBBBBB 
;BLBBBBIBBB 
SBOBBBBBBBB 
“BB lg BBIBBB 

Sbebbbbbbbb 
SBEBBBBBBBB 
*BQBBBBBBBB 
EBI3BC3BBBBBB 
;iiniimi 
iBEBBBBBBflfl 
iBKXBBg 

saisssl I 

uiibSSII 

BBISBBBIBBEBBI 

BBKciBC!BG9B[£BU3l 

fiBIEBKt)B[!3Bl3BI2l 

B— BEBBIBIHI 


IBBISBEBBI 

■■■■■■■■I 

IBBEBBBBI 


lEBEBElBEIBEiBEBEII 

lEBIIMniBBBn 

lICBEBISBEaBeniSBEII 

IISBBBBBBBBBBBBI 
lEflflflfllBBBIB fll 


ICZBlHI 

IBIQI 


ICSI 


IBEHHMI 

IBEBISIkSinSIBlIBBBI 
IBEBBBEBBBBBBBBI 
■EBBBBBBBBBBBBI 

IBB 

IStBB 

■E BBBBBBBBBBBBI 


::nsHC!BEHdBiE:ac 3 aBi 

!■&■■■■■■■■■■■■! 
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Li sti ng of No. 12 Data Cards for the DUGLFT Test Problem 


Columns : 

11 

21 

32 34 




i 

1 

i 1 





0.5 

2 

WSQF 

7 


0.5 

0.5 


WSOF 

8 


0.5 



WSQF 

9 

0.5 


0.5 

1 

WSQF 

10 

0.5 

0.5 

0.5 


WSQF 

1 1 

0.5 

0.5 



WSQF 

12 

1 . 0 


0.5 

1 

WSQF 

13 

1 .0 

0.5 

0.5 


WSQF 

14 

1 . 0 

0.5 



WSQF 

15 

1 . 5 


0.5 

1 

WSQF 

16 

1 .5 

0.5 

0.5 


WSQF 

17 

1 .5 

0.5 

0.0427 


WSQF 

18 

2.0 


0.5 

1 

WSQF 

19 

2.0 

0.5 

0.5 


WSQF 

20 

2.0 

0.5 



WSQF 

21 

2.5 


0.5 

1 

WSQF 

22 

2.5 

0.5 

0.5 


WSQF 

23 

2 . 5 

0.5 



WSQF 

24 

3.0 


0.5 

1 

WSQF 

25 

3.0 

0.5 

0.5 


WSQF 

26 

3.0 

0.5 



WSQF 

27 


0.5 


2 

WSQF 

28 


0.5 

-0.5 


WSQF 

29 



-0.5 


WSQF 

30 

0.5 

0.5 


1 

WSQF 

31 

0.5 

0.5 

-0.5 


WSQF 

32 

0.5 


-0.5 


WSQF 

33 

1 .0 

0.5 


1 

WSQF 

34 

1 .0 

0.5 

-0.5 


WSQF 

35 

1 .0 


-0.5 


WSQF 

36 

1 .5 

0.5 

-0.0427 

1 

WSQF 

37 

1 .5 

0.5 

-0.5 


WSQF 

38 

1 .5 


-0.5 


WSQF 

39 

2.0 

0.5 


1 

WSQF 

40 

2.0 

0.5 

-0.5 


WSQF 

41 

2.0 


-0.5 


WSQF 

42 

2.5 

0.5 


1 

WSQF 

43 

2.5 

0.5 

-0.5 


WSQF 

44 

2.5 


-0.5 


WSQF 

45 

3.0 

0.5 


1 

WSQF 

46 

3.0 

0.5 

-0.5 


WSQF 

47 

3 . 0 


— 0 . 5 


WSQF 

48 

2. 

3.5 


2 1 

WSQF 

49 

1 .85 

3.5 

-.0134 


WSQF 

50 

1 .6 

3.5 

-.0353 


WSQF 

51 

1 .3 

3.5 

-.05 


WSQF 

52 

1 .1 

3.5 

-.0361 


WSQF 

53 

1 . 

3.5 



WSQF 

54 

1 .1 

3.5 

.0361 


WSQF 

55 

1 .3 

3.5 

.05 


WSQF 

56 

1 .6 

3.5 

.0353 


WSQF 

57 

1 .85 

3.5 

.0134 


WSQF 

58 

2.0 

3.5 



WSQF 

59 

2.15 

3.5 



WSQF 

60 

2.3 

3.5 



WSQF 

61 

2.5 

3.5 



WSQF 

62 

2. 

2.5 


1 

WSQF 

63 

1 .85 

2.5 

-.0134 


WSQF 

64 

1 .6 

2.5 

-.0353 


WSQF 

65 

1 .3 

2.5 

-.05 


WSQF 

66 


110 



1 . t 

2.5 

-.0361 


WSQF 

67 

1 .0 

2.5 



WSQF 

68 

1 .1 

2 . 5 

.0361 


WSQF 

69 

1 .3 

2.5 

.05 


WSQF 

70 

1 .6 

2.5 

.0353 


WSQF 

71 

1 .85 

2.5 

.0134 


WSQF 

72 

2. 

2 . 5 



WSQF 

73 

2.15 

2.5 



WSQF 

74 

2.3 

2.5 



WSQF 

75 

2.5 

2.5 



WSQF 

76 

2.0 

1 .5 


1 

WSQF 

77 

1 .85 

1 .5 

-.0134 


WSQF 

78 

1 .6 

1 .5 

-.0353 


WSQF 

79 

1 .3 

1 . 5 

-.05 


WSQF 

60 

1.1 

1 .5 

-.0361 


WSQF 

81 

1.0 

1 .5 



WSQF 

82 

1.1 

1.5 

.0361 


WSQF 

83 

1 .3 

1 .5 

. 05 


WSQF 

84 

1 .6 

1 .5 

.0353 


WSQF 

85 

1 .85 

1 .5 

.0134 


WSQF 

86 

2.0 

1 .5 



WSQF 

87 

2.15 

1 .5 



WSQF 

88 

2.3 

1 .5 



WSQF 

89 

2.5 

1.5 



WSQF 

90 

2.0 

.5 


1 

WSQF 

91 

1.85 

.5 

-.0134 


WSQF 

92 

1.6 

.5 

-.0353 


WSQF 

93 

1.3 

.5 

-.05 


WSQF 

94 

1 .1 

.5 

-.0361 


WSQF 

95 

1 .0 

.5 



WSQF 

96 

1 . 1 

.5 

.0361 


WSQF 

97 

1.3 

.5 

.05 


WSQF 

98 

1 .6 

.5 

.0353 


WSQF 

99 

1.85 

.5 

.0134 


WSQF 

100 

2.0 

. 5 



WSQF 

101 

2.15 

.5 



WSQF 

102 

2.3 

.5 



WSQF 

103 

2.5 

.5 



WSQF 

104 

2.0 



1 

WSQF 

105 

1 .85 


-.0134 


WSQF 

106 

1 .6 


-.0353 


WSQF 

107 

1 .3 


-.05 


WSQF 

108 

1.1 


-.0361 


WSQF 

109 

1.0 




WSQF 

110 

1.1 


.0361 


WSQF 

111 

1.3 


.05 


WSQF 

112 

1 .6 


• 0353 


WSQF 

113 

1.85 


.0134 


WSQF 

114 

2.0 




WSQF 

115 

2.15 




WSQF 

116 

2.3 




WSQF 

117 

2.5 



3 

WSOF 

118 


in 
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o 

00 

■1 
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■i 

si 

■I 

w 



■1 
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3 
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El 


HI 
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THE FORCE COMPONENTS OF THE SECTION ARE -.241962E+00 0. .274433E+01 

THE MOMENT COMPONENTS OF THE SECTION ARE .518401E+01 -.331023E+01 .437360E+00 
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THE MOMENT COMPONENTS OF THE ENTIRE BODY ARE .516046E+01 -.364205E+01 .563370E+00 


PROGRAM DUGLFT DOUGLAS AIRCRAFT COMPANY PAGE 25. 

MODIFIED BY H.G.NORMENT, ATMOSPHERIC SCIENCE ASSOCIATES* BEDFORD* MASS. 
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END OF THE VELOCY ROUTINE 
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INDICATES THE POINT IS INSIDE THE BODY 



APPENDIX A 


PRANDTL-GLAUERT COMPRESSIBILITY CORRECTION 


Sears (ref. 15) shows that application of the small perturbation 
approximation to Euler's equation for subsonic flow results in 


2 

3 3 > 


xx 


+ $ + $ 

yy 


zz 


where the velocity potential function $ is given by 


(Al) 


* = “Ux - <J> (A2) 

U is uniform free stream speed in the positive x axis direction, «j> is per- 
turbation velocity potential function, and if N M is free stream Mach number, 
then 


» 2 ■ 1 - (A3) 

Subscripts indicate partial differentiation, for example, $ xx = 9 2 $/9x9x. 
In terms of perturbation velocity (u, v, w), eq. (Al ) is 

e 2 u x + v y + w z = 0 . (A4) 

If the body surface is defined by the function S(x,y,z) = 0, then the 
surface boundary condition is satisfied by 

(U + u)S x + vS y + wS z = 0. (A5) 

The Prandtl-Glauert transformation stretches the coordinate system 
(x,y,z) and velocity (u,v,w) to {x,y, z) and (u,v,w) according to 

x = x /3 

y = y (A6) 

z = z 
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and 


2 

a = g u 

v = gv (A7) 

w = gw 

and the transformed partial derivatives of S are 
5 * = pS x 

S = S (A8) 

y y 


In terms of the transformed quantities eqs. (A4) and (A5) are 

u + v + w = 0 (A9) 

x y z 


and 


(U + u/g )S^ + vSy + miS z - 0 (A10) 

Use and significance of these results in the context of the Hess flow 
calculation method are as follows. We apply the eq. (A6) transformation to 
all body surface coordinates as they are first input to the code. This 
amounts to stretching the body along the x axis by a factor of 1/g. Veloci- 
ties (u, v , w) then are calculated for the stretched body, and source 
strengths are determined for the quadrilateral panels subject to the 
boundary condition of eq. (A10). From eq. (A9) we see that these veloci- 
ties satisfy the continuity equation for potential flow. Finally, when 
velocities in the original coordinate system are recovered, via application 
of the inverse transformation 
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(All) 


2 

u = a/e 
v = v/e 
w = w/e 

they are already corrected as they stand for significant compressibility 
effects caused by the high Mach number free stream flow. 

It should be understood that this compressibility correction is based on 
application of the small perturbation approximation. As such it assumes 
that perturbation velocity, pressure and density, and their derivatives, 
are small relative to the free stream. Thus, for example, the method does 
not treat compressibility effects caused by locally high Mach number flow. 

In addition to the transformation above, there are others that are 
needed to calculate pressure force and moments on the body. Pressure 
force, F, on a quadrilateral element is given by 

F = - nC p A (A12) 

where n is the unit vector normal to the element surface, A is the element 
area, and 



->- ->■ 

where V is flow velocity at the element centroid. Both n and A must be 
o 

calculated for the original, unstretched body, and they are recovered from 
quantities computed for the stretched body as follows. 

Change our subscript notation such that (n x , n^, n^) represents com- 
ponents of n in the x, y and z axis directions. Then we have 

"x ' "x /c 

n * n B/c 
j y 

n z = n z e/c 
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where 


2 2 2 2 2 

c = ft + 0 (n + n) 
x y z 


( A1 5 ) 


Using the notation of Appendix B below, area for nonlifting elements, 
A^, and for lifting elements A^, are calculated by the equations 


^NL - ^ 


1 ' ' 2 


no - n, 


)/2 


( A1 6 ) 


and 


A^ ( ?2 ” £3 ~ 5^)(n^ - 0^)72 


(A17) 


These equations can be evaluated from the following general equations for 
the £ and n differences 


n - h = 5 


Si ( n ■ x j> + a i2 { »i - i'j 1 + a i3 (z i - z j> 


(A18) 


n i ' n j ^ 


B a 21 ( Xi - Xj) + (B n z a u - - y.) 


+ { \ a 12 - H )(z i " Z J ] 


(A19) 


where 


2 _ 2 2 , 2 2 

q B a ll + a 12 + a 13 


(A20) 


and the are are transformation matrix elements for transforming coordinates 
between the flow and element reference frames (see ref. 4, sec. 9.2). 
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APPENDIX B 


BODY PENETRATION TESTS 


INTRODUCTION 

A sequence of three sets of tests are applied to each on-body surface 
element for which the exact source and dipole contribution equations are 
required to determine if a particle trajectory has penetrated the element 
surface. In the event that the tests of the first two sequences (see 
pp. 67-68) fail to eliminate the possibility of penetration, then this 
means that the plane of the element has been intersected by the trajectory 
during the current integration time step, and it is the purpose of the 
sequence 3 tests, which are described in detail in this appendix, to de- 
termine if the intersection point lies within or on the boundaries of the 
element. 

Since the particle is sufficiently close to the element centroid that 
the "exact" equations are being used for the velocity contribution calcula- 
tions (see ref. 3, sec. 7.4 and ref. 4, sec. 9.52), all coordinates have 
been transformed to the element reference frame: origin at the element 

centroid, x and y axes in the plane of the element with the z axis normal 
to the element plane and positive toward the exterior of the body. Element 
corner points are defined differently in the element coordinate systems for 
nonlifting and lifting elements, and therefore, we treat these types of 
elements separately below. 

First we find the line that connects the current particle position 
with its position at the previous time step, and find the point of inter- 
section of this line in the element plane. Designate the current particle 
coordinates as (x c , y c , z c ) and the coordinates at the previous time step 
as (x r , y , z r ). Then the line that connects these points may be defined 
by the ratios. 
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X - 


X 


r 




(Bl) 


and since in the element coordinate system the element plane is defined 
by z = 0, the coordinates of the point of intersection, (x , y , z ), of 
this line with the element plane are ^ 


x = 
P 


z (x - x )/(z - z ) 

c r c r c 


y p = y c - z c <*r • *c )/(z r - 2 c> 


z = 0 
P 


NONLIFTING ELEMENTS 


(B2) 


Nonlifting elements are configured as shown in Figure Bl in which the 
corner points are numbered in accordance with the convention shown on p. 30 
above. The element coordinate system is defined such that points 1 and 3 
lie on the x axis; point 2 or 4, but not both, also may lie on the x axis. 

To determine if point P (Fig. B2) is within or on the element boundaries, 
consider the m and n*" neighboring element corner points, and define vectors 
v mn between these points and v m p between the points m and P as 


v 


mn 






(B3) 


v mP 



(B4) 


where i and j are unit vectors in the directions of the x and y axes, and 
the pair (m,n) has values (1,2), (2,3), (3,4) or (4,1). If we take the 
vector product of v mn with k, the unit vector in the direction of the z 
axis, we obtain a vector, ffr mn , in the element plane that is normal to v mn 
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X 


Figure Bl. Configuration of a nonlifting element in 
the element coordinate system. The posi- 
tive z axis direction is up from the paper. 



I 

I 


Figure B2. Relationship of element corner points m and 
n, the point P of intersection of the par- 
ticle trajectory with the element plane, and 
vectors v„„, v„ n and N . 


and is directed toward the interior of the element, given by 



<”n - - (E n 



(B5) 


Now, if the scalar product, S mn , of N mn with v m p is negative, the point P 
must be exterior to the element. Thus we compute 


S mn = < X P - ‘ Tn> ' % - %,)<«„ - « m ) < B6 > 


successively for the (m,n) pairs given above, and if the relation 


S 


mn 


< 0 


(B7) 


is valid for any (m,n), point P cannot lie on or within the element bounda- 
ries. If relation (B7) fails for all of the (m,n), we have proved that 
P lies on or within the element boundaries, and thus, we have proved that 
the trajectory intersects the body surface. 


LIFTING ELEMENTS 

Lifting elements are configured in the element coordinate system as 
shown in Figure B3, with sides 1-2 and 3-4 parallel with the x axis. Thus, 
it is immediately apparent that point P cannot lie inside or on the 
element boundaries if (y - oi ) > 0 or (y - n-J < 0. These tests are 

pi p O 

combined in the code such that if 

(y p - n-, ) (y p - n 3 ) > 0 (B8) 

there can be no impaction. 
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Figure B3. Configuration of a lifting element in the element 
coordinate system. The positive z axis is directed 
up from the paper. The x axis is approximately 
parallel with the streamwise direction, and the 
y axis is approximately parallel with the spanwise 
direction of an airfoil. 


If relation (B8) is not satisfied, P lies on or between the line 
defined by points 1 and 2 and the line defined by points 3 and 4. Then 
we proceed to see which sides of the lines 2-3 and 4-1 P lies on by 
applying essentially the same method as was used for the nonlifting 
element. In terms of the slope, M = ax/ a y , of lines 2-3 and 4-1, we 
test the relations 


M 23 <*p - *3> ' X p + S 3 •= 0 

(B9) 

M 4i % - h > - * P - > 0 

(BIO) 


and if either is satisfied, P cannot lie inside or on the element bounda- 
ries. If neither is satisifed, we have proved that P lies inside or on 
the boundaries of the element, and thus, we have proved that the trajec- 
tory has intersected the body surface. 
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APPENDIX C 


SPECIAL SOURCE AND DIPOLE CONTRIBUTION CALCULATIONS 
INTRODUCTION 

Unless precautions are taken numerical problems can arise when velocities 
are calculated at points very close to element surfaces, to element edges, 
and for lifting elements, at points close to extensions of element edges. 

This is because some of the equation terms for exact calculation of source 
and dipole contributions become infinite at such points. In this appendix 
we identify the troublesome terms, and describe how the numerical problems 
are handled in the velocity calculations. The notations of Hess (ref. 3) 
and Hess and Smith (ref. 4) are used. Our frame of reference is the element 
coordinate system in which nonlifting and lifting elements are defined as 
illustrated in Figures B1 and B3 above. Additional definitions are given in 
Figure Cl. 

GENERAL 

If the coordinates of the point for which the velocity calculation is 
desired are (x, y, z), then the distance of this point from the i^* 1 element 
corner point is 


r . = ^(x - + (y - o^ 2 + z 2 (Cl) 

and for certain purposes the calculation point coordinates relative to the 
corner point are defined as 

«, = (X - 5,1/r, 

= (y - n i )/r i (C2) 

Y, - z /**i 
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If r. is very small equations (C2) may blow up, therefore if r. is calculated 
by equation (Cl) to be less than 10 -3 t, where t is the maximum diagonal of 
the element, then we set r. = 10 3 t. That is, we impose the restriction 


(r • ) = 10" 3 t. 

min 


(C3) 


H 


(ij) 

k 


Define quantity h£ 1J ), which applies to nonlifting (NL) 
which applies to lifting (L) elements, as 


elements, and 


h 


(ij) 

k 


m ij a k " 6 k 


(C4NL) 


where 


m 


ij 



(C5NL) 


and 


H 


(ij) 

k 


= M. .g. 
ij k 


a k 


(C4L) 


where 



(C5L) 


and k - i or k = j. For equations (C4L) and (C5L) we are restricted to 
(ij) pairs (32) and (41). (See Figures B1 and B3.) 

It is easy to show that, as the projection of point (x, y, z) in the 
z = 0 plane approaches the straight line that passes through both the i^ 
and j th element corner points, h^'^for ti[^) approaches zero. For lifting 
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elements, as the projection of point (x, y, z) in the z = 0 plane approaches 
one of the straight lines that pass through points 1 and 2 or 3 and 4, we 
have -* 0. 

Finally, it is obvious that as point (x, y, z) approaches the element 
plane, z 0, and usually -* 0. 

SOURCE CONTRIBUTION TERMS 

The Hess-Smith method of velocity calculation requires that potential 
sources be uniformly distributed over the surface of each on-body element. 

In calculating source velocity contributions from individual elements 
numerical problems are encountered at points very close to or on element 
edges, and at general points close to or in the plane of the element. 

The x component of velocity contributed by the distributed source 
on an element consists of a summation of terms of the form 

Hi ~ n i L ( i j ) 

d. . 
ij 

where 


L (ij) 



+ 



) 


and the y component consists of a sum of terms of the form 



5. 




L (ij) 


where 


(C6) 



V 


<v 


(,) 2 + (l 


-" 1 >‘ 


(C7) 
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is the distance between the neighboring i^ and element corner points. 
From equation (C6) it is apparent that as point (x, y, z) approaches the 
line connecting the i th and j th points, Thus, we impose the 


restriction that the minimum value of (r. + r. 

3 i j 

10 : that is 


d ij )/(r i + r j + d ij) be 


(L (ij) ) . = Jin(10" 3 ) 

min 

(C8) 

Contribution to the z component of velocity from an 
summation of terms 

element consists 


(C9NL) 

= tan-‘(Vk 2 + 6 k H k 1J> \ 

(C9L) 


Thus as the point (x, y, z) approaches the plane of the element, Yk -* 0 
ai j|d -*• + w/2. If a k h^ 1J ^ (or also approaches zero, then 

Tr J is indeterminate, but the summation is such that indeterminate 
terms of equal magnitude but opposite sign cancel, so that we simply set 
the T^) to zero. In summary, the following restrictions are imposed. 


T (ij) 

'k 




<W 1j) 


> e 


T (ij) 

'k 


= 0 


irt\, < £ > 


a h ( i J ) 
k n k 

o H (ij ) 
e k H k 


< e 


(C10a) 


(ClOb) 
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where e is specified to be 10~ 4 for nonlifting elements and 10" 3 for lifting 
elements: the different values being imposed more for convenience than 

necessity. In equation (ClOa) the sign of tt/2 is taken to be the sign of 
the product Y k a k h^ 1J ^ (or Y k e k H k i: ^)- 

It can be shown that the sum of terms results in the following 

values of the z component of velocity, v , contributed by the element for 
the following calculation points which are a small distance e from the 
plane of the element: 


Location of (x, y, e) 

1. Inside the boundaries of the element 

2. On a line connecting two element cor- 
ner points 

3. On an element corner point 

4. Outside the element boundaries 


(sign of e) 2ir 
(sign of e) * 

(sign of e) tt/2 
0 


The restrictions imposed by equations (ClOa) and (Cl Ob) automatically 
ensure that these results are obtained. 


DIPOLE CONTRIBUTION TERMS 

Hess uses distributed dipoles to represent lift vorticity sheets, and 
he finds that space partial derivatives of various moments of the dipole 
distributions are required in computing velocity contributions from the 
lift vorticity. Numerical problems may arise in computing space deriva- 
tives of L^^, which are required for the evaluation of the space deriva- 
tives of the and <t> 10 potential dipole moments, and in computing space 
derivatives of T^ J '^ , which are required to evaluate space derivatives of 
the 4 >qq potential dipole moment. 

The derivatives of L^ J ^ are 


3L (1 J> ...... 3 L (1 J> 

— — D . a . + a . ) , — - 

3x 1J 1 3 3y 


= D..(e. + $.), — - — = D . . (v . + y • ) 
ij i j az ij VT i Tr 


(ij) 
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where 


/ 


D u ■ 




(r i+ rj) 2 -d 2 


(Cll) 


J ' 

As point (x, y, z) approaches the line connecting element corner points i 
and j , D. . blov 

' J 

striction that 


and j, D.. blows up. Consistent with equation (C8) we impose the re- 
l J 


( r i + r j^ " d ij 2 = 10_3 ^ r ’ + r < + d -s-i) 2 • 


min 


i J U 


(ij) 


Space derivatives of J; are given by 


3T, 


(ij) 


n I e k + “k p k 


(ij) 


8X 




where 


3 T k (, ' j) 

- Yk ( 

' 6 kki - 

p k (iJ) ] 

+ H k (1 i> 

3y 

■ nr \ 

k ^ 


2 

aT 

9 k 

- — ( 

[m.. p 
L ij - 

k (,J) ] ■ 

\ " B k H k 

az 

' r k \ 

. ^ 


2 

p (ij) _ 
K k 

M ij Y k + 

R H ^ ij ’^ 
B k H k 




(ij)' 


(Cl 2) 


(C13a) 


(Cl 3b) 


(C13c) 


(C14) 
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It is easily seen for points (x, y, e), where e is small or zero, that are 
not on or very close to element boundary lines or their extensions, (i.e.,, 

| 3 k | »0 and »0), that equations (C13) can be used as given without 

any trouble. If the calculation point approaches an extension of an element 
edge line (i.e., y k 0, 3 k -> 0 or H^ 1J ^ + 0, and r? + r? - d?.» 0), then 
special equations can be used. However, if an element edge is approached, 
the only practical response seems to be to set the partial derivatives of 
eq. (C13) to zero, as is shown below. 

The general equation for the derivative of <}>qq with respect to space 
coordinates ft (= x, y or z) is 



The terms in this equation are most conveniently considered in pairs as 
is indicated by the For example. 



and so forth. In the following we assume that the line that passes through 
quadrilateral points (ij) is very closely approached, where (ij) is the super- 
script on ip. The quantities 6 and e are as defined in Figure Cl. 

Case 1. Calculation point approaches an extension of an element edge: 

^(ij) _ q f or n _ x anc | y. 

^ 12) = ( M 41 I x ‘ 5 1 I _1 ‘ M 32 I x ’ I "0 Sin20 
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^ (^32 | x ~^3 | ” ^41 | x “ ^4 | j ® 

♦ ™ = M ij ( |y ■ n j -1 - | y - n i | _1 j Sin 1 2 e; (ij) = (32) or (41) 

(Note: The factor sin 2 e is missing from these equations, (7.8.5), in refer- 
ence 3. ) 

Case 2. Calculation point approaches an element edge: 

Jij) = o for (ij) = (12) and (34). 

A 

6 appears in the denominator of the equations for all other 

To accommodate these problems without seriously compromising overall 
accuracy, for example for particle trajectory calculations, we have selected 
the following sequence of tests which are applied in order: 

1. If fz/t | >_ 10~ 3 , bypass the following tests and use equations 

( Cl 3 ) . Here, as elsewhere, t is the maximum diagonal of the 
element. 

2. If |r k 6 k /t | >_ 10” 3 and |r k H k ^ 1J Vtj >_ 10~ 3 use equations (C13). 

3. If both of the above tests are negative, then set 
9T (ij) 3T( ij ’) 3T 

K K K _ r\ 
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One additional special calculation is required for semi -infinite 
wake elements (see section 7.10 of reference 3) which may be needed in 
calculating space derivatives of such as 


g L ( 34 ) _ a 4 - 1 
9X r 4 - (x - e 4 ) 


(C16) 


To prevent blow-up when the denominator becomes very small, we impose the 
restriction 


[ 





t 


(C17 


where t is the maximum diagonal of the element. 
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